0.1 Preprocessing:

Basic dada2 pipeline for 16S:
truncLen=c(200, 180)
trimLeft=c(19,20)
maxEE=c(3,4)
truncQ=2
pool=pseudo
database - SILVA 138.1

Delete all bad taxa with parameters below:
annosighted on the phylum level +
chloroplasts/mithondria +
low reads/richness(based on the plot results)

tree from SEPP-QIIME2-plagin:

qiime tools import
–input-path ref_filt.fasta
–output-path rep_seq.qza
–type ‘FeatureData[Sequence]’ qiime fragment-insertion sepp
–i-representative-sequences rep_seq.qza   –i-reference-database sepp-refs-silva-128.qza
–o-tree insertion-tree.qza
–o-placements insertion-placements.qza
–p-threads 40
unzip -p insertion-tree.qza /data/ > tree.nwk

ITS preprocessing:

dada2 with cutadapt trimming
maxEE=c(5,8)
truncQ = 8
minLen = 50
pool=“pseudo”
UNITE version - 10.05.2021

ML tree from IQ-TREE programm:

mafft –auto FASTA > FASTA_ALIGHNED
iqtree -s FASTA_ALIGHNED -nt 30 -B 1000

0.2 Import

Import phyloseq object and libraries
Dataset - straw decomposition time series - from factor Day - 10 levels Also in the dataset added bulk soils samples.

library(phyloseq)
library(tidyverse)
library(ggpubr)
library(ampvis2)
library(ANCOMBC)
library(heatmaply)
library(compositions)
library(igraph)
library(WGCNA)
require(DESeq2)
require(phyloseq)

ps_bulk <- readRDS("ps_bulk")

0.3 functions

#phyloseq object to ampvis2 object
#https://gist.github.com/KasperSkytte/8d0ca4206a66be7ff6d76fc4ab8e66c6

phyloseq_to_ampvis2 <- function(physeq) {
  #check object for class
  if(!any(class(physeq) %in% "phyloseq"))
    stop("physeq object must be of class \"phyloseq\"", call. = FALSE)
  
  #ampvis2 requires taxonomy and abundance table, phyloseq checks for the latter
  if(is.null(physeq@tax_table))
    stop("No taxonomy found in the phyloseq object and is required for ampvis2", call. = FALSE)
  
  #OTUs must be in rows, not columns
  if(phyloseq::taxa_are_rows(physeq))
    abund <- as.data.frame(phyloseq::otu_table(physeq)@.Data)
  else
    abund <- as.data.frame(t(phyloseq::otu_table(physeq)@.Data))
  
  #tax_table is assumed to have OTUs in rows too
  tax <- phyloseq::tax_table(physeq)@.Data
  
  #merge by rownames (OTUs)
  otutable <- merge(
    abund,
    tax,
    by = 0,
    all.x = TRUE,
    all.y = FALSE,
    sort = FALSE
  )
  colnames(otutable)[1] <- "OTU"
  
  #extract sample_data (metadata)
  if(!is.null(physeq@sam_data)) {
    metadata <- data.frame(
      phyloseq::sample_data(physeq),
      row.names = phyloseq::sample_names(physeq), 
      stringsAsFactors = FALSE, 
      check.names = FALSE
    )
    
    #check if any columns match exactly with rownames
    #if none matched assume row names are sample identifiers
    samplesCol <- unlist(lapply(metadata, function(x) {
      identical(x, rownames(metadata))}))
    
    if(any(samplesCol)) {
      #error if a column matched and it's not the first
      if(!samplesCol[[1]])
        stop("Sample ID's must be in the first column in the sample metadata, please reorder", call. = FALSE)
    } else {
      #assume rownames are sample identifiers, merge at the end with name "SampleID"
      if(any(colnames(metadata) %in% "SampleID"))
        stop("A column in the sample metadata is already named \"SampleID\" but does not seem to contain sample ID's", call. = FALSE)
      metadata$SampleID <- rownames(metadata)
      
      #reorder columns so SampleID is the first
      metadata <- metadata[, c(which(colnames(metadata) %in% "SampleID"), 1:(ncol(metadata)-1L)), drop = FALSE]
    }
  } else
    metadata <- NULL
  
  #extract phylogenetic tree, assumed to be of class "phylo"
  if(!is.null(physeq@phy_tree)) {
    tree <- phyloseq::phy_tree(physeq)
  } else
    tree <- NULL
  
  #extract OTU DNA sequences, assumed to be of class "XStringSet"
  if(!is.null(physeq@refseq)) {
    #convert XStringSet to DNAbin using a temporary file (easiest)
    fastaTempFile <- tempfile(pattern = "ampvis2_", fileext = ".fa")
    Biostrings::writeXStringSet(physeq@refseq, filepath = fastaTempFile)
  } else
    fastaTempFile <- NULL
  
  #load as normally with amp_load
  ampvis2::amp_load(
    otutable = otutable,
    metadata = metadata,
    tree = tree,
    fasta = fastaTempFile
  )
}

#variance stabilisation from DESeq2

ps_vst <- function(ps, factor){
  

  
  diagdds = phyloseq_to_deseq2(ps, as.formula(paste( "~", factor)))              
  diagdds = estimateSizeFactors(diagdds, type="poscounts")
  diagdds = estimateDispersions(diagdds, fitType = "local") 
  pst <- varianceStabilizingTransformation(diagdds)
  pst.dimmed <- t(as.matrix(assay(pst))) 
  # pst.dimmed[pst.dimmed < 0.0] <- 0.0
  ps.varstab <- ps
  otu_table(ps.varstab) <- otu_table(pst.dimmed, taxa_are_rows = FALSE) 
  return(ps.varstab)
} 

#WGCNA visualisation
#result - list class object with attributes:
# ps - phyloseq object
# amp - ampwis2 object
# heat - heatmap with absalute read numbers
# heat_rel - hetmap with relative abundances
# tree - phylogenetic tree with taxonomy

color_filt <- function(ps, df){
  
      library(tidyverse)
      library(reshape2)
      library(gridExtra)

  
  l = list()
  for (i in levels(df$module)){
    message(i)
    color_name <-  df %>% 
      filter(module == i) %>% 
      pull(asv) %>% 
      unique()
    ps.col <- prune_taxa(color_name, ps)
    amp.col <- phyloseq_to_ampvis2(ps.col)
    heat <- amp_heatmap(amp.col, tax_show = 60, 
              group_by = "Day", 
              tax_aggregate = "OTU",
              tax_add = "Genus", 
              normalise=FALSE, 
              showRemainingTaxa = TRUE)
    
    ps.rel  <-  phyloseq::transform_sample_counts(ps, function(x) x / sum(x) * 100)
    ps.rel.col <- prune_taxa(color_name, ps.rel)
    amp.r <- phyloseq_to_ampvis2(ps.rel.col)
    heat.rel <- amp_heatmap(amp.r, tax_show = 60, 
          group_by = "Day", 
          tax_aggregate = "OTU",
          tax_add = "Genus", 
          normalise=FALSE, 
          showRemainingTaxa = TRUE)
    
    tree <- ps.col@phy_tree
    taxa <- as.data.frame(ps.col@tax_table@.Data)
    p1 <- ggtree(tree) + 
              geom_tiplab(size=2, align=TRUE, linesize=.5) +
              theme_tree2()
    
    taxa[taxa == "Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium"] <- "Allorhizobium"
    taxa[taxa == "Burkholderia-Caballeronia-Paraburkholderia"] <- "Burkholderia"
  
    tx <- taxa %>% 
      rownames_to_column("id") %>% 
      mutate(id = factor(id, levels = rev(get_taxa_name(p1)))) %>% 
      dplyr::select(-c(Kingdom, Species, Order)) %>% 
      melt(id.var = 'id')

    p2 <- ggplot(tx, aes(variable, id)) + 
      geom_tile(aes(fill = value), alpha = 0.4) +
      geom_text(aes(label = value), size = 3) +
      theme_bw() + 
      theme(legend.position = "none",
             axis.ticks.x=element_blank(),
             axis.text.y=element_blank(),
             axis.ticks.y=element_blank(),
             axis.title.x = element_blank(),
            axis.title.y = element_blank()) 

    p <- ggpubr::ggarrange(p1, p2, widths = c(0.9, 1))
    l[[i]] <- list("ps" = ps.col, 
                   "amp" = amp.col,
                   "heat" = heat,
                   "heat_rel" = heat.rel,
                   "tree" = p,
                   "taxa" = knitr::kable(taxa))
                    
  }
  return(l)
}


color_filt_broken <- function(ps, df, ps.pruned){
  
      library(tidyverse)
      library(reshape2)
      library(gridExtra)

  
  l = list()
  for (i in levels(df$module)){
    message(i)
    color_name <-  df %>% 
      filter(module == i) %>% 
      pull(asv) %>% 
      unique()
    ps.col <- prune_taxa(color_name, ps)
    ps.col.pruned <- prune_taxa(color_name, ps.pruned)
    amp.col <- phyloseq_to_ampvis2(ps.col)
    heat <- amp_heatmap(amp.col, tax_show = 60, 
              group_by = "Day", 
              tax_aggregate = "OTU",
              tax_add = "Genus", 
              normalise=FALSE, 
              showRemainingTaxa = TRUE)
    
    ps.rel  <-  phyloseq::transform_sample_counts(ps.col, function(x) x / sum(x) * 100)
    amp.r <- phyloseq_to_ampvis2(ps.rel)
    heat.rel <- amp_heatmap(amp.r, tax_show = 60, 
          group_by = "Day", 
          tax_aggregate = "OTU",
          tax_add = "Genus", 
          normalise=FALSE, 
          showRemainingTaxa = TRUE)
    
    tree <- ps.col.pruned@phy_tree
    taxa <- as.data.frame(ps.col.pruned@tax_table@.Data)
    p1 <- ggtree(tree) + 
              geom_tiplab(size=2, align=TRUE, linesize=.5) +
              theme_tree2()
    
    taxa[taxa == "Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium"] <- "Allorhizobium"
    taxa[taxa == "Burkholderia-Caballeronia-Paraburkholderia"] <- "Burkholderia"
  
    tx <- taxa %>% 
      rownames_to_column("id") %>% 
      mutate(id = factor(id, levels = rev(get_taxa_name(p1)))) %>% 
      dplyr::select(-c(Kingdom, Species, Order)) %>% 
      melt(id.var = 'id')

    p2 <- ggplot(tx, aes(variable, id)) + 
      geom_tile(aes(fill = value), alpha = 0.4) +
      geom_text(aes(label = value), size = 3) +
      theme_bw() + 
      theme(legend.position = "none",
             axis.ticks.x=element_blank(),
             axis.text.y=element_blank(),
             axis.ticks.y=element_blank(),
             axis.title.x = element_blank(),
            axis.title.y = element_blank()) 

    p <- ggpubr::ggarrange(p1, p2, widths = c(0.9, 1))
    l[[i]] <- list("ps" = ps.col, 
                   "amp" = amp.col,
                   "heat" = heat,
                   "heat_rel" = heat.rel,
                   "tree" = p,
                   "taxa" = knitr::kable(taxa))
                    
  }
  return(l)
}

detachAllPackages <- function() {

  basic.packages <- c("package:stats","package:graphics","package:grDevices","package:utils","package:datasets","package:methods","package:base")

  package.list <- search()[ifelse(unlist(gregexpr("package:",search()))==1,TRUE,FALSE)]

  package.list <- setdiff(package.list,basic.packages)

  if (length(package.list)>0)  for (package in package.list) detach(package, character.only=TRUE, force = TRUE)

}

plot_alpha_w_toc <- function(ps, group, metric) {
  
  require(phyloseq)
  require(ggplot2)
  
  ps_a <- prune_taxa(taxa_sums(ps) > 0, ps)
  
  er <- estimate_richness(ps_a)
  df_er <- cbind(ps_a@sam_data, er)
  df_er <- df_er %>% select(c(group, metric))
  stat.test <- aov(as.formula(paste0(metric, "~", group)), data = df_er) %>%
    rstatix::tukey_hsd()
  y <-  seq(max(er[[metric]]), length=length(stat.test$p.adj), by=max(er[[metric]]/20))

  plot_richness(ps_a, x=group, measures=metric) + 
    geom_boxplot() +
    geom_point(size=1.2, alpha=0.3) +
    ggpubr::stat_pvalue_manual(
      stat.test, 
      label = "p.adj.signif", 
      y.position = y) +
    theme_light() + 
    scale_color_brewer(palette="Dark2") +
    theme(axis.text.x = element_text(angle = 90),
          axis.title.x=element_blank()) +
    labs(y=paste(metric, "index")) 
}


# standart NMDS plot tool frop phyloseq with some additional aestatics 
# have stress value on plot - may work as fuck
beta_custom_norm_NMDS_elli_w <- function(ps, seed = 7888, normtype="vst", Color="What", Group="Repeat"){
  require(phyloseq)
  require(ggplot2)
  require(ggpubr)
  library(ggforce)
  
  
  ordination.b <- ordinate(ps, "NMDS", "bray")
  mds <- as.data.frame(ordination.b$points)
  p  <-  plot_ordination(ps,
                         ordination.b,
                         type="sample",
                         color = Color,
                         title="NMDS - Bray-Curtis",
                         # title=NULL,
                         axes = c(1,2) ) + 
    theme_bw() + 
    theme(text = element_text(size = 10)) + 
    geom_point(size = 3) +
    annotate("text",
    x=min(mds$MDS1) + abs(min(mds$MDS1))/7,
    y=max(mds$MDS2),
    label=paste0("Stress -- ", round(ordination.b$stress, 3))) +
    geom_mark_ellipse(aes_string(group = Group, label = Group),
                      label.fontsize = 10,
                      label.buffer = unit(2, "mm"),
                      label.minwidth = unit(5, "mm"),
                      con.cap = unit(0.1, "mm"),
                      con.colour='gray') +
    theme(legend.position = "none") +
    ggplot2::scale_fill_viridis_c(option = "H")
  
  return(p)
}


# alpha with aov + tukie post-hock - useless, but it looks pretty good
plot_alpha_w_toc <- function(ps, group, metric) {
  
  require(phyloseq)
  require(ggplot2)
  
  ps_a <- prune_taxa(taxa_sums(ps) > 0, ps)
  
  er <- estimate_richness(ps_a)
  df_er <- cbind(ps_a@sam_data, er)
  df_er <- df_er %>% select(c(group, metric))
  stat.test <- aov(as.formula(paste0(metric, "~", group)), data = df_er) %>%
    rstatix::tukey_hsd()
  y <- seq(max(er[[metric]]), length=length(stat.test$p.adj.signif[stat.test$p.adj.signif != "ns"]), by=max(er[[metric]]/20))
  
  plot_richness(ps_a, x=group, measures=metric) + 
    geom_boxplot() +
    geom_point(size=1.2, alpha=0.3) +
    stat_pvalue_manual(
      stat.test, 
      label = "p.adj.signif", 
      y.position = y,
      hide.ns=TRUE) +
    theme_light() + 
    scale_color_brewer(palette="Dark2") +
    theme(axis.text.x = element_text(angle = 90),
          axis.title.x=element_blank()) +
    labs(y=paste(metric, "index")) 
}

lsf.str()
## beta_custom_norm_NMDS_elli_w : function (ps, seed = 7888, normtype = "vst", Color = "What", Group = "Repeat")  
## color_filt : function (ps, df)  
## color_filt_broken : function (ps, df, ps.pruned)  
## detachAllPackages : function ()  
## phyloseq_to_ampvis2 : function (physeq)  
## plot_alpha_w_toc : function (ps, group, metric)  
## ps_vst : function (ps, factor)

0.3.1 modify metadata


Add Group parameter to metadata -
- early - D01, D03, D05
- middle - D07, D08, D10
- late - D13, D14, D15

sample.data <- ps_bulk@sam_data %>% 
  data.frame() %>% 
  mutate(Group = if_else(Day %in% c("D01", "D03", "D05"), "early", 
                         if_else(Day %in% c("D07", "D08","D10"), "middle", 
                                 if_else(Day %in% c("D12", "D13", "D14", "D15"), "late", "bulk")))) %>% 
  dplyr::rename('Bag' = 'Day') %>% 
  mutate(Group = factor(Group, levels=c("early", "middle","late", "bulk"))) %>% 
  mutate(Day = Bag %>% 
           forcats::fct_recode( "3" = "D01",
                                "14" = "D03",
                                "28" = "D05",
                                "49" = "D07" ,
                                "63" = "D08",
                                "91" = "D10",
                                "119" = "D12",
                                "140" = "D13",
                                "161" = "D14",
                                "182" = "D15")
         ) %>% 
  phyloseq::sample_data()

sample_data(ps_bulk) <- sample.data

 ps_bulk@sam_data %>% 
  DT::datatable(caption = "Sample data of phyloseq - 16S")

0.3.2 Select from the dataset only experemental part

ps.f <- prune_samples(!sample_data(ps_bulk)$Day %in% "bulk", ps_bulk)
ps.f  <- prune_taxa(taxa_sums(ps.f) > 0, ps.f)
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
ps.f@sam_data  %>% 
  DT::datatable(caption = "Sample data of phyloseq - 16S")

0.4 Alpha diversity

With bulk soil

p1 <- plot_alpha_w_toc(ps = ps_bulk, group = "Day", metric = "Observed")
## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
##   # Was:
##   data %>% select(group)
## 
##   # Now:
##   data %>% select(all_of(group))
## 
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
##   # Was:
##   data %>% select(metric)
## 
##   # Now:
##   data %>% select(all_of(metric))
## 
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
p2 <- plot_alpha_w_toc(ps = ps_bulk, group = "Day", metric = "Shannon")
p3 <- plot_alpha_w_toc(ps = ps_bulk, group = "Day", metric = "InvSimpson")
ggarrange(p1, p2, p3, nrow = 1)

Without bulk soil

p1 <- plot_alpha_w_toc(ps = ps.f, group = "Day", metric = "Observed")
p2 <- plot_alpha_w_toc(ps = ps.f, group = "Day", metric = "Shannon")
p3 <- plot_alpha_w_toc(ps = ps.f, group = "Day", metric = "InvSimpson")
ggarrange(p1, p2, p3, nrow = 1)

p.observed <- plot_alpha_w_toc(ps = ps.f, group = "Group", metric = c("Observed")) + 
  theme(axis.title.y = element_blank())
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
p.shannon <- plot_alpha_w_toc(ps = ps.f, group = "Group", metric = c("Shannon")) + 
  theme(axis.title.y = element_blank())
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
p.simpson <- plot_alpha_w_toc(ps = ps.f, group = "Group", metric = c("InvSimpson")) + 
  theme(axis.title.y = element_blank())
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
ggpubr::ggarrange(p.observed, p.shannon, p.simpson, ncol = 3)

0.4.0.1 mpd

mpd - index of alpha diversity “tree shagginess” - counted separately because it is in a separate package with built-in validity (based on permutations)
here is a link about the meanin of values in the columns:
https://www.rdocumentation.org/packages/picante/versions/1.8.2/topics/ses.mpd
in general the early stage is significantly less diverse than the other stages

physeq_merged <- merge_samples(ps.f, "Group", fun=sum)
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
# ps.f@sam_data

# picante::mpd(l_vst$blue$ps@otu_table@.Data, cophenetic(l_vst$blue$ps@phy_tree)) %>% 
#   mean(na.rm = TRUE)
# 
# picante::mpd(l_vst$salmon$ps@otu_table@.Data, cophenetic(l_vst$salmon$ps@phy_tree)) %>% 
#   mean(na.rm = TRUE)

mpd.res <- picante::ses.mpd(physeq_merged@otu_table@.Data, cophenetic(physeq_merged@phy_tree)) 

as.data.frame(mpd.res)
##        ntaxa  mpd.obs mpd.rand.mean mpd.rand.sd mpd.obs.rank  mpd.obs.z
## early    344 1.707396      1.840030  0.02736870            1 -4.8462083
## middle   445 1.832365      1.840085  0.02325692          377 -0.3319760
## late     687 1.852563      1.841628  0.01467574          769  0.7451306
##        mpd.obs.p runs
## early      0.001  999
## middle     0.377  999
## late       0.769  999

0.5 beta statistic

0.5.1 permanova

permanova - group significantlly different - dispersion between is more, than inside groups

dist <- phyloseq::distance(ps.f, "bray")
metadata <- as(sample_data(ps.f@sam_data), "data.frame")
vegan::adonis2(dist ~ Group, data = metadata)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = dist ~ Group, data = metadata)
##          Df SumOfSqs      R2      F Pr(>F)    
## Group     2   3.2510 0.33952 8.2247  0.001 ***
## Residual 32   6.3243 0.66048                  
## Total    34   9.5753 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

0.5.2 dynamic of beta diversity

Bray Curtis distance by steps between each day.   (For example, all D15 versus D14.) Can we also put soil respiration on this picture as well?

ps.f.r <- rarefy_even_depth(ps.f, rngseed = 777)

avg.r <- ps.f.r@otu_table %>% 
  as.data.frame() %>%
  vegan::avgdist(10)

avg <- ps.f@otu_table %>% 
  as.data.frame() %>%
  vegan::avgdist(10)

# avg %>%
#   as.matrix() %>%
#   as_tibble(rownames= "sample") %>%
#   pivot_longer(-sample) %>%
#   filter(sample < name) %>%
#   mutate(repeat_a = str_replace(sample, ".*-", ""),
#          repeat_b = str_replace(name, ".*-", ""), 
#          day_a = as.numeric(str_replace(sapply(strsplit(sample, "-"), `[`, 3), "D", "")),
#           day_b = as.numeric(str_replace(sapply(strsplit(name, "-"), `[`, 3), "D", "")),
#          diff = abs(day_a - day_b),
#          early = day_a < 10) %>% 
#   filter(repeat_a == repeat_b & diff < 10) %>%
#   group_by(diff, repeat_a, early) %>%
#   summarize(median = median(value)) %>%
#   ungroup() %>%
#   ggplot(aes(x=diff, y=median, color=early, group=paste0(repeat_a, early))) +
#   geom_line(size=0.25) +
#   geom_smooth(aes(group=early), se=FALSE, size=4) +
#   labs(x="Distance between time points",
#        y="Median Bray-Curtis distance") +
#   scale_x_continuous(breaks=1:9) +
#   scale_color_manual(name=NULL,
#                      breaks=c(TRUE, FALSE),
#                      values=c("blue", "red"),
#                      labels=c("Early", "Late")) +
#   guides(color = guide_legend(override.aes = list(size=1))) +
#   theme_classic()

avg.r %>%
  as.matrix() %>%
  as_tibble(rownames= "sample") %>%
  pivot_longer(-sample) %>%
  filter(sample < name) %>%
  mutate(repeat_a = str_replace(sample, ".*-", ""),
         repeat_b = str_replace(name, ".*-", ""), 
         day_a = as.numeric(str_replace(sapply(strsplit(sample, "\\."), `[`, 3), "D", "")),
         day_b = as.numeric(str_replace(sapply(strsplit(name, "\\."), `[`, 3), "D", ""))) %>% 
  mutate(day_a = as.numeric(as.factor(day_a) %>% forcats::fct_recode("2" = "3", "3" = "5", "4" = "7", "5" = "8", "6" = "10", "7" = "13", "8" = "14", "9" = "14", "10" = "15")),
         day_b = as.numeric(as.factor(day_b) %>% forcats::fct_recode("2" = "3", "3" = "5", "4" = "7", "5" = "8", "6" = "10", "7" = "12", "8" = "13", "9" = "14", "10" = "15")),
         diff = abs(day_a - day_b)) %>% 
  filter(diff == 1) %>% 
  mutate(day_b = as.factor(day_b) %>% forcats::fct_recode("14" = "2", "28" = "3","49" = "4",  "63" = "5", "91" = "6", "161" = "9", "140" = "8", "182" = "10", "119" = "7")) %>% 
  ggplot(aes(x=day_b, y=value)) +
    geom_boxplot() +
  theme_bw() +
  labs(x="Days",
       y="Bray-Curtis distance")

0.5.3 basic NMDS plot

From the picture it looks like there is a gradual slowing down of changes -
but from the previous picture it follows that this is not the case - there is a dip between 7-8-10,
but on average the points differ quite consistently.

beta_custom_norm_NMDS_elli_w(ps.f, C="Group", G="Day")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1063746 
## Run 1 stress 0.1063746 
## ... New best solution
## ... Procrustes: rmse 0.00002091149  max resid 0.00007447643 
## ... Similar to previous best
## Run 2 stress 0.1427868 
## Run 3 stress 0.1099872 
## Run 4 stress 0.1396556 
## Run 5 stress 0.1063746 
## ... Procrustes: rmse 0.00000856638  max resid 0.00004570378 
## ... Similar to previous best
## Run 6 stress 0.1063746 
## ... Procrustes: rmse 0.00001753402  max resid 0.00006418496 
## ... Similar to previous best
## Run 7 stress 0.1411972 
## Run 8 stress 0.1063746 
## ... Procrustes: rmse 0.00004854709  max resid 0.0001716062 
## ... Similar to previous best
## Run 9 stress 0.1099872 
## Run 10 stress 0.1063499 
## ... New best solution
## ... Procrustes: rmse 0.00155448  max resid 0.006797445 
## ... Similar to previous best
## Run 11 stress 0.1448551 
## Run 12 stress 0.1099721 
## Run 13 stress 0.1482809 
## Run 14 stress 0.1438721 
## Run 15 stress 0.1099872 
## Run 16 stress 0.1409441 
## Run 17 stress 0.1393472 
## Run 18 stress 0.1099873 
## Run 19 stress 0.143872 
## Run 20 stress 0.1099872 
## *** Best solution repeated 1 times
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation ideoms with `aes()`

0.6 Splitting the dataset

0.6.1 Basic information about dataset

amp <- phyloseq_to_ampvis2(ps.f)
amp
## ampvis2 object with 5 elements. 
## Summary of OTU table:
## [1]       35     1063   510548     4412    32228    12752 14587.09
## 
## Assigned taxonomy:
##      Kingdom       Phylum        Class        Order       Family        Genus 
##   1063(100%)   1063(100%)   1063(100%)   1063(100%) 1054(99.15%)  736(69.24%) 
##      Species 
##    90(8.47%) 
## 
## Metadata variables: 5 
##  SampleID, Bag, Description, Group, Day
amp_heatmap(amp,
            tax_show = 22,
            group_by = "Day",
            tax_aggregate = "Phylum",
            tax_class = "Proteobacteriota",
            normalise=TRUE,
            showRemainingTaxa = TRUE)
## Warning: Transformation introduced infinite values in discrete y-axis

amp_heatmap(amp,
            tax_show = 40,
            group_by = "Day",
            tax_aggregate = "OTU",
            tax_add = "Genus",
            normalise=TRUE,
            showRemainingTaxa = TRUE)

Let’s divide the dataset into two groups.

1st - more than 10% of samples should contain at least 10 reads - this group will be analyzed further

ps.inall <- phyloseq::filter_taxa(ps.f, function(x) sum(x > 10) > (0.1*length(x)), TRUE)
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
amp.inall <- phyloseq_to_ampvis2(ps.inall)
amp_heatmap(amp.inall,
            tax_show = 40,
            group_by = "Day",
            tax_aggregate = "OTU",
            tax_add = "Genus",
            normalise=FALSE,
            showRemainingTaxa = TRUE)

The second part - the remaining, but more than 100 reeds in all samples(hereinafter - outliers majors(OM))
The remaining phylotypes are dropped from the analysis.

ps.exl  <- phyloseq::filter_taxa(ps.f, function(x) sum(x > 10) < (0.1*length(x)), TRUE)
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
ps.exl  <- prune_taxa(taxa_sums(ps.exl) > 100, ps.exl)
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
amp.exl <- phyloseq_to_ampvis2(ps.exl)
amp_heatmap(amp.exl,
            tax_show = 40,
            group_by = "Day",
            tax_aggregate = "OTU",
            tax_add = "Genus",
            normalise=FALSE,
            showRemainingTaxa = TRUE)

Same, but with relative abundance.

ps.per  <-  phyloseq::transform_sample_counts(ps.f, function(x) x / sum(x) * 100) 
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
ps.exl.taxa <- taxa_names(ps.exl)
ps.per.exl <- prune_taxa(ps.exl.taxa, ps.per)
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
amp.exl.r <- phyloseq_to_ampvis2(ps.per.exl)

amp_heatmap(amp.exl.r, tax_show = 60, 
      group_by = "Day", 
      tax_aggregate = "OTU",
      tax_add = "Genus", 
      round =  2, 
      normalise=FALSE, 
      showRemainingTaxa = TRUE)

amp_heatmap(amp.exl.r, tax_show = 60, 
      group_by = "Day", 
      tax_aggregate = "OTU",
      tax_add = "Genus", 
      round =  2, 
      normalise=FALSE, )

OM aggregated by genus. Relative abundance.

amp_heatmap(amp.exl.r,
            tax_show = 30, 
            group_by = "Day", 
            tax_aggregate = "Genus",
            tax_add = "Phylum",
            tax_class = "Proteobacteria",
            round =  2,
            normalise=FALSE, 
            showRemainingTaxa = TRUE)

If we look at how OMs are distributed by day - it seems that the distribution of OMs is related not only to the specificity of individual bags, but also to purely technical features - the outlier values are mostly not with biological repeats, but with technical ones (red line - selected visually)

p_box <- phyloseq::sample_sums(ps.per.exl) %>% 
  as.data.frame(col.names = "values") %>% 
  setNames(., nm = "values") %>% 
  rownames_to_column("samples") %>% 
  mutate(Day = sapply(strsplit(samples, "\\."), `[`, 3)) %>% 
  ggplot(aes(x=Day, y=values, color=Day, fill = Day)) +
  geom_boxplot(aes(color=Day, fill = Day)) +
  geom_point(color = "black", position = position_dodge(width=0.2)) +
  geom_hline(yintercept = 10, colour = "red") +
  theme_bw() +
  theme(legend.position = "none")
     
p_box <- p_box + viridis::scale_color_viridis(option = "H", discrete = TRUE, direction=1, begin=0.1, end = 0.9, alpha = 0.5)

p_box + viridis::scale_fill_viridis(option = "H", discrete = TRUE, direction=1, begin=0.1, end = 0.9, alpha = 0.3)

0.7 WGCNA

0.7.1 clr

after clr normalization / looks disgusting - let’s try vst normalization from DESeq2

otu.inall <- as.data.frame(ps.inall@otu_table@.Data) 
ps.inall.clr <- ps.inall
otu_table(ps.inall.clr) <- phyloseq::otu_table(compositions::clr(otu.inall), taxa_are_rows = FALSE)
data <- ps.inall.clr@otu_table@.Data %>% 
  as.data.frame()
rownames(data) <- as.character(ps.inall.clr@sam_data$Description)
powers <-  c(c(1:10), seq(from = 12, to=30, by=1))
sft <-  pickSoftThreshold(data, powerVector = powers, verbose = 5, networkType = "signed")
## pickSoftThreshold: will use block size 321.
##  pickSoftThreshold: calculating connectivity for given powers...
##    ..working on genes 1 through 321 of 321
## Warning: executing %dopar% sequentially: no parallel backend registered
##    Power SFT.R.sq  slope truncated.R.sq  mean.k. median.k.  max.k.
## 1      1    0.133  19.60         0.8930 160.0000 160.00000 164.000
## 2      2    0.309 -18.70         0.8910  83.1000  82.90000  90.300
## 3      3    0.465 -11.70         0.8350  44.9000  44.50000  52.900
## 4      4    0.362  -5.99         0.7610  25.2000  24.60000  32.600
## 5      5    0.252  -3.34         0.7890  14.6000  14.00000  21.000
## 6      6    0.218  -2.19         0.7400   8.7500   8.31000  14.100
## 7      7    0.257  -1.82         0.7160   5.4300   5.09000   9.800
## 8      8    0.328  -1.71         0.6800   3.4900   3.21000   7.030
## 9      9    0.497  -1.92         0.7830   2.3100   2.09000   5.380
## 10    10    0.629  -2.03         0.7780   1.5800   1.39000   4.270
## 11    12    0.231  -4.27         0.0333   0.7990   0.66100   2.900
## 12    13    0.861  -2.22         0.9010   0.5900   0.47100   2.450
## 13    14    0.850  -2.16         0.8560   0.4450   0.34100   2.110
## 14    15    0.262  -3.70         0.0795   0.3430   0.24700   1.840
## 15    16    0.265  -3.55         0.0859   0.2690   0.18400   1.620
## 16    17    0.267  -3.42         0.0912   0.2140   0.13700   1.430
## 17    18    0.902  -1.97         0.9060   0.1730   0.10500   1.280
## 18    19    0.912  -1.92         0.9100   0.1420   0.08000   1.150
## 19    20    0.924  -1.87         0.9200   0.1170   0.06170   1.040
## 20    21    0.292  -4.09         0.1340   0.0981   0.04810   0.944
## 21    22    0.304  -3.06         0.1950   0.0829   0.03820   0.861
## 22    23    0.300  -2.94         0.1830   0.0706   0.03050   0.788
## 23    24    0.294  -2.79         0.1680   0.0607   0.02460   0.736
## 24    25    0.928  -1.62         0.9080   0.0525   0.01950   0.693
## 25    26    0.935  -1.60         0.9170   0.0458   0.01530   0.656
## 26    27    0.933  -1.57         0.9140   0.0402   0.01210   0.623
## 27    28    0.264  -2.43         0.0682   0.0354   0.00967   0.593
## 28    29    0.266  -2.40         0.0708   0.0314   0.00769   0.565
## 29    30    0.196  -2.67        -0.0255   0.0280   0.00608   0.541
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n", main = paste("Scale independence"))
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], labels=powers,cex=0.9,col="red")
abline(h=0.9,col="salmon")

0.7.2 vst

after vst normalisation

ps.varstab <- ps_vst(ps.inall, "Day")

data2 <- ps.varstab@otu_table@.Data %>% 
  as.data.frame()

rownames(data2) <- as.character(ps.varstab@sam_data$Description)
powers <-  c(seq(from = 1, to=10, by=0.5), seq(from = 11, to=20, by=1))
sft2 <-  pickSoftThreshold(data2, powerVector = powers, verbose = 5, networkType = "signed hybrid")
## pickSoftThreshold: will use block size 321.
##  pickSoftThreshold: calculating connectivity for given powers...
##    ..working on genes 1 through 321 of 321
##    Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
## 1    1.0    0.339 -0.985         0.7120  46.800  43.10000  90.50
## 2    1.5    0.569 -1.150         0.8650  28.800  25.60000  65.50
## 3    2.0    0.711 -1.200         0.9340  18.900  16.90000  49.10
## 4    2.5    0.754 -1.210         0.9080  13.000  11.20000  38.10
## 5    3.0    0.807 -1.280         0.9260   9.260   7.36000  30.40
## 6    3.5    0.864 -1.330         0.9610   6.820   5.08000  24.80
## 7    4.0    0.863 -1.350         0.9310   5.160   3.59000  20.60
## 8    4.5    0.845 -1.380         0.8890   3.990   2.56000  17.30
## 9    5.0    0.780 -1.450         0.8240   3.140   1.84000  14.70
## 10   5.5    0.814 -1.440         0.8690   2.520   1.35000  12.70
## 11   6.0    0.819 -1.440         0.8580   2.050   1.03000  11.00
## 12   6.5    0.889 -1.350         0.9420   1.690   0.78200   9.66
## 13   7.0    0.868 -1.340         0.9070   1.410   0.60600   8.53
## 14   7.5    0.873 -1.330         0.9020   1.190   0.46700   7.68
## 15   8.0    0.818 -1.360         0.8180   1.020   0.38800   6.97
## 16   8.5    0.833 -1.360         0.8250   0.877   0.32000   6.35
## 17   9.0    0.854 -1.360         0.8550   0.762   0.25600   5.83
## 18   9.5    0.892 -1.370         0.9150   0.666   0.21400   5.37
## 19  10.0    0.886 -1.350         0.8940   0.587   0.17900   4.96
## 20  11.0    0.958 -1.290         0.9790   0.465   0.11800   4.29
## 21  12.0    0.902 -1.290         0.9060   0.376   0.08010   3.75
## 22  13.0    0.922 -1.270         0.9270   0.310   0.05370   3.31
## 23  14.0    0.913 -1.240         0.9080   0.260   0.03730   2.94
## 24  15.0    0.936 -1.190         0.9260   0.221   0.02670   2.63
## 25  16.0    0.924 -1.180         0.9110   0.190   0.02050   2.36
## 26  17.0    0.181 -1.840        -0.0518   0.165   0.01440   2.16
## 27  18.0    0.910 -1.210         0.8860   0.145   0.00980   2.05
## 28  19.0    0.901 -1.210         0.8770   0.129   0.00664   1.95
## 29  20.0    0.909 -1.240         0.8970   0.115   0.00468   1.87
plot(sft2$fitIndices[,1], -sign(sft2$fitIndices[,3])*sft2$fitIndices[,2], xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n", main = paste("Scale independence"))
text(sft2$fitIndices[,1], -sign(sft2$fitIndices[,3])*sft2$fitIndices[,2], labels=powers,cex=0.9,col="red")
abline(h=0.9,col="salmon")

# WGCNA, as old and not supported
detachAllPackages()
library(WGCNA)


net3 <- WGCNA::blockwiseModules(data2,
                          power=5.5,
                          TOMType="signed",
                          networkType="signed hybrid",
                          nThreads=0)


mergedColors2 <-  WGCNA::labels2colors(net3$colors, colorSeq = c("salmon", "darkgreen", "cyan", "red", "blue", "plum"))

plotDendroAndColors(
  net3$dendrograms[[1]],
  mergedColors2[net3$blockGenes[[1]]],
  "Module colors",
  dendroLabels = FALSE,
  hang = 0.03,
  addGuide = TRUE,
  guideHang = 0.05)

unique(ps.inall@sam_data$Bag)
##  [1] D01 D03 D05 D07 D08 D10 D12 D13 D14 D15
## Levels: D01 D03 D05 D07 D08 D10 D12 D13 D14 D15

0.7.2.1 import libraries 2

library(phyloseq)
library(tidyverse)
library(ggpubr)
library(ampvis2)
library(heatmaply)
library(WGCNA)
library(phyloseq)
library(ggtree)
library(tidyverse)
library(KneeArrower)
modules_of_interest = mergedColors2 %>% 
  unique()

module_df <- data.frame(
  asv = names(net3$colors),
  colors = mergedColors2
)

# module_df[module_df == "yellow"] <- "salmon"

submod <-  module_df %>%
  subset(colors %in% modules_of_interest)

row.names(module_df) = module_df$asv

subexpr = as.data.frame(t(data2))[submod$asv,]

submod_df <- data.frame(subexpr) %>%
  mutate(
    asv = row.names(.)
  ) %>%
  pivot_longer(-asv) %>%
  mutate(
    module = module_df[asv,]$colors
  )

submod_df <- submod_df %>% 
  mutate(name = gsub("\\_.*","",submod_df$name)) %>% 
  group_by(name, asv) %>% 
  summarise(value = mean(value), asv = asv, module = module) %>% 
  relocate(c(asv, name, value, module)) %>% 
  ungroup() %>% 
  mutate(module = as.factor(module))
## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
##   always returns an ungrouped data frame and adjust accordingly.
p <- submod_df %>% 
  ggplot(., aes(x=name, y=value, group=asv)) +
  geom_line(aes(color = module),
            alpha = 0.2) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 90),
    legend.position = "none") +
  facet_grid(rows = vars(module)) +
  labs(x = "day",
       y = "normalized abundance")

p + scale_color_manual(values = levels(submod_df$module)) +
  scale_x_discrete(labels=(c("D01"="3",
                                "D03"="14",
                                "D05"="28",
                                "D07"="49" ,
                                "D08"="63",
                                "D10"="91",
                                "D12"="119",
                                "D13"="140",
                                "D14"="161",
                                "D15"="182")
         ))

0.7.2.2 ordination

UNIFRAC - the sausage narrows down

ps.inall.col <- ps.inall
df <- module_df %>% 
  rename("id" = "asv")
df <- df %>% 
  dplyr::select(-"id") %>% 
  mutate(colors = as.factor(colors))
taxa <- as.data.frame(ps.inall@tax_table@.Data)
tx <- cbind(taxa, df)
tx$colors <- factor(tx$colors, levels = c("salmon", "darkgreen", "cyan", "red", "blue", "plum"))
tax_table(ps.inall.col) <- tax_table(as.matrix(tx))
ord <- ordinate(ps.inall.col, "NMDS", "unifrac")
## Run 0 stress 0.07519118 
## Run 1 stress 0.05788475 
## ... New best solution
## ... Procrustes: rmse 0.04301752  max resid 0.1702092 
## Run 2 stress 0.07049655 
## Run 3 stress 0.07832129 
## Run 4 stress 0.07112217 
## Run 5 stress 0.08037211 
## Run 6 stress 0.05788483 
## ... Procrustes: rmse 0.00006064862  max resid 0.0002316986 
## ... Similar to previous best
## Run 7 stress 0.07519005 
## Run 8 stress 0.07078907 
## Run 9 stress 0.07050821 
## Run 10 stress 0.05777393 
## ... New best solution
## ... Procrustes: rmse 0.005923683  max resid 0.02753314 
## Run 11 stress 0.08098516 
## Run 12 stress 0.0577739 
## ... New best solution
## ... Procrustes: rmse 0.00004120185  max resid 0.0001533557 
## ... Similar to previous best
## Run 13 stress 0.07358919 
## Run 14 stress 0.05777396 
## ... Procrustes: rmse 0.0002527585  max resid 0.0009829406 
## ... Similar to previous best
## Run 15 stress 0.0577738 
## ... New best solution
## ... Procrustes: rmse 0.0000951983  max resid 0.0003642492 
## ... Similar to previous best
## Run 16 stress 0.05788472 
## ... Procrustes: rmse 0.005926721  max resid 0.02756207 
## Run 17 stress 0.08029901 
## Run 18 stress 0.07056578 
## Run 19 stress 0.08054233 
## Run 20 stress 0.07358869 
## *** Best solution repeated 1 times
plot_ordination(ps.inall.col, ord, type = "species", color = "colors") + 
  scale_color_manual(values =  c("salmon", "darkgreen", "cyan", "red", "blue", "plum")) +
  theme_bw() +
  theme(legend.position = "none")

Bray - the sausage even
Late stages on the left - with early clusters overlapping, late clusters seems separated What affects on the axis2? There is clearly some sort of pattern.

ord <- ordinate(ps.inall.col, "NMDS", "bray")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.09266802 
## Run 1 stress 0.09266802 
## ... Procrustes: rmse 0.000007178359  max resid 0.00001887827 
## ... Similar to previous best
## Run 2 stress 0.09266802 
## ... Procrustes: rmse 0.000008206866  max resid 0.00002707994 
## ... Similar to previous best
## Run 3 stress 0.09266802 
## ... Procrustes: rmse 0.00001652581  max resid 0.00006336731 
## ... Similar to previous best
## Run 4 stress 0.09266802 
## ... Procrustes: rmse 0.000009393617  max resid 0.00003872337 
## ... Similar to previous best
## Run 5 stress 0.09266804 
## ... Procrustes: rmse 0.00004581879  max resid 0.0001753545 
## ... Similar to previous best
## Run 6 stress 0.09266803 
## ... Procrustes: rmse 0.00003403899  max resid 0.0001493913 
## ... Similar to previous best
## Run 7 stress 0.09266802 
## ... New best solution
## ... Procrustes: rmse 0.000002174393  max resid 0.00001011168 
## ... Similar to previous best
## Run 8 stress 0.1298078 
## Run 9 stress 0.1387607 
## Run 10 stress 0.1383668 
## Run 11 stress 0.09266803 
## ... Procrustes: rmse 0.00001333261  max resid 0.00004836488 
## ... Similar to previous best
## Run 12 stress 0.09266805 
## ... Procrustes: rmse 0.00004695051  max resid 0.0001789074 
## ... Similar to previous best
## Run 13 stress 0.1297207 
## Run 14 stress 0.09266802 
## ... Procrustes: rmse 0.000003347503  max resid 0.00000917534 
## ... Similar to previous best
## Run 15 stress 0.1296264 
## Run 16 stress 0.09266802 
## ... Procrustes: rmse 0.000004806095  max resid 0.00001233449 
## ... Similar to previous best
## Run 17 stress 0.1296263 
## Run 18 stress 0.09266803 
## ... Procrustes: rmse 0.00001123112  max resid 0.00004466975 
## ... Similar to previous best
## Run 19 stress 0.1298078 
## Run 20 stress 0.09266803 
## ... Procrustes: rmse 0.00001845328  max resid 0.00006284917 
## ... Similar to previous best
## *** Best solution repeated 7 times
plot_ordination(ps.inall.col, ord, type = "species", color = "colors") + 
  scale_color_manual(values =  c("salmon", "darkgreen", "cyan", "red", "blue", "plum")) +
  theme_bw() +
  theme(legend.position = "none")

0.7.2.3 visualisation

Next are the same pictures for all the groups.

l_vst <- color_filt(ps.inall, submod_df)
l_vst

$cyan \(cyan\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 82 taxa and 35 samples ] sample_data() Sample Data: [ 35 samples by 4 sample variables ] tax_table() Taxonomy Table: [ 82 taxa by 7 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 82 tips and 81 internal nodes ] refseq() DNAStringSet: [ 82 reference sequences ]

\(cyan\)amp ampvis2 object with 5 elements. Summary of OTU table: [1] 35 82 70590 183 5339 1720 2016.86

Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 82(100%) 82(100%) 82(100%) 82(100%) 82(100%) 71(86.59%) 12(14.63%)

Metadata variables: 5 SampleID, Bag, Description, Group, Day

\(cyan\)heat \(cyan\)heat_rel \(cyan\)tree \(cyan\)taxa

Kingdom Phylum Class Order Family Genus Species
Seq143 Bacteria Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae Burkholderia NA
Seq65 Bacteria Pseudomonadota Gammaproteobacteria Burkholderiales Alcaligenaceae NA NA
Seq149 Bacteria Pseudomonadota Gammaproteobacteria Burkholderiales Comamonadaceae Variovorax NA
Seq503 Bacteria Cyanobacteriota Vampirivibrionia Obscuribacterales Obscuribacteraceae NA NA
Seq622 Bacteria Pseudomonadota Alphaproteobacteria Ferrovibrionales Ferrovibrionaceae Ferrovibrio soli
Seq341 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Xanthobacteraceae NA NA
Seq239 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Beijerinckiaceae Methylobacterium-Methylorubrum NA
Seq232 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Beijerinckiaceae Bosea thiooxidans
Seq538 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Beijerinckiaceae Bosea NA
Seq105 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Xanthobacteraceae Tardiphaga robiniae
Seq31 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Xanthobacteraceae Starkeya NA
Seq212 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Hyphomicrobiaceae Hyphomicrobium NA
Seq605 Bacteria Pseudomonadota Alphaproteobacteria Micropepsales Micropepsaceae NA NA
Seq386 Bacteria Pseudomonadota Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingopyxis NA
Seq610 Bacteria Pseudomonadota Alphaproteobacteria Sphingomonadales Sphingomonadaceae Altererythrobacter NA
Seq595 Bacteria Pseudomonadota Alphaproteobacteria Sphingomonadales Sphingomonadaceae Altererythrobacter NA
Seq132 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Rhizobiaceae Mesorhizobium NA
Seq462 Bacteria Pseudomonadota Alphaproteobacteria Caulobacterales Caulobacteraceae NA NA
Seq77 Bacteria Actinobacteriota Thermoleophilia Solirubrobacterales Solirubrobacteraceae Conexibacter NA
Seq549 Bacteria Myxococcota Polyangia Polyangiales Polyangiaceae Pajaroellobacter NA
Seq345 Bacteria Myxococcota Polyangia Polyangiales Polyangiaceae Pajaroellobacter NA
Seq103 Bacteria Myxococcota Polyangia Polyangiales Polyangiaceae Sorangium NA
Seq486 Bacteria Actinobacteriota Acidimicrobiia Microtrichales Ilumatobacteraceae NA NA
Seq113 Bacteria Actinobacteriota Actinobacteria Micrococcales Promicromonosporaceae Promicromonospora NA
Seq598 Bacteria Actinobacteriota Actinobacteria Micrococcales Microbacteriaceae Leifsonia NA
Seq171 Bacteria Actinobacteriota Actinobacteria Micrococcales Microbacteriaceae NA aoyamense
Seq674 Bacteria Actinobacteriota Actinobacteria Micrococcales Microbacteriaceae Agromyces NA
Seq139 Bacteria Actinobacteriota Actinobacteria Corynebacteriales Mycobacteriaceae Mycobacterium NA
Seq159 Bacteria Actinobacteriota Actinobacteria Streptosporangiales Streptosporangiaceae Herbidospora mongoliensis
Seq654 Bacteria Actinobacteriota Actinobacteria Streptosporangiales Streptosporangiaceae Nonomuraea NA
Seq86 Bacteria Actinobacteriota Actinobacteria Streptosporangiales Thermomonosporaceae Actinocorallia NA
Seq267 Bacteria Actinobacteriota Actinobacteria Streptomycetales Streptomycetaceae Streptomyces NA
Seq745 Bacteria Actinobacteriota Actinobacteria Propionibacteriales Propionibacteriaceae Jiangella NA
Seq245 Bacteria Bacteroidota Bacteroidia Sphingobacteriales Sphingobacteriaceae Pedobacter panaciterrae
Seq189 Bacteria Bacteroidota Bacteroidia Sphingobacteriales Sphingobacteriaceae Pedobacter NA
Seq492 Bacteria Bacteroidota Bacteroidia Sphingobacteriales NS11-12 marine group NA NA
Seq716 Bacteria Bacteroidota Bacteroidia Sphingobacteriales NS11-12 marine group NA NA
Seq309 Bacteria Bacillota Bacilli Bacillales Bacillaceae Bacillus NA
Seq101 Bacteria Bacillota Bacilli Paenibacillales Paenibacillaceae Paenibacillus NA
Seq52 Bacteria Bacillota Bacilli Paenibacillales Paenibacillaceae Paenibacillus NA
Seq18 Bacteria Bacillota Bacilli Paenibacillales Paenibacillaceae Paenibacillus NA
Seq265 Bacteria Bacillota Bacilli Paenibacillales Paenibacillaceae Paenibacillus NA
Seq415 Bacteria Bacillota Bacilli Paenibacillales Paenibacillaceae Paenibacillus NA
Seq439 Bacteria Bacillota Bacilli Bacillales Planococcaceae Domibacillus NA
Seq121 Bacteria Bacillota Bacilli Bacillales Bacillaceae Bacillus NA
Seq54 Bacteria Bacillota Bacilli Bacillales Planococcaceae Lysinibacillus NA
Seq144 Bacteria Bacillota Bacilli Bacillales Planococcaceae Lysinibacillus NA
Seq84 Bacteria Verrucomicrobiota Verrucomicrobiae Chthoniobacterales Terrimicrobiaceae Terrimicrobium NA
Seq349 Bacteria Verrucomicrobiota Verrucomicrobiae Chthoniobacterales Terrimicrobiaceae Terrimicrobium NA
Seq228 Bacteria Verrucomicrobiota Verrucomicrobiae Chthoniobacterales Terrimicrobiaceae Terrimicrobium NA
Seq631 Bacteria Verrucomicrobiota Verrucomicrobiae Chthoniobacterales Chthoniobacteraceae Chthoniobacter NA
Seq451 Bacteria Verrucomicrobiota Verrucomicrobiae Chthoniobacterales Chthoniobacteraceae Chthoniobacter NA
Seq431 Bacteria Verrucomicrobiota Verrucomicrobiae Verrucomicrobiales Verrucomicrobiaceae Roseimicrobium gellanilyticum
Seq315 Bacteria Verrucomicrobiota Verrucomicrobiae Verrucomicrobiales Verrucomicrobiaceae Verrucomicrobium spinosum
Seq246 Bacteria Planctomycetota Planctomycetes Gemmatales Gemmataceae Gemmata NA
Seq375 Bacteria Planctomycetota Planctomycetes Isosphaerales Isosphaeraceae Singulisphaera NA
Seq395 Bacteria Acidobacteriota Acidobacteriae Acidobacteriales Acidobacteriaceae (Subgroup 1) Edaphobacter NA
Seq73 Bacteria Bacteroidota Bacteroidia Cytophagales Microscillaceae Ohtaekwangia NA
Seq237 Bacteria Bacteroidota Bacteroidia Cytophagales Microscillaceae Ohtaekwangia NA
Seq420 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Pseudoflavitalea NA
Seq32 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Pseudoflavitalea NA
Seq765 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Niastella NA
Seq40 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Niastella hibisci
Seq234 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Niastella NA
Seq372 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Niastella NA
Seq126 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Flavitalea NA
Seq155 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Chitinophaga ginsengihumi
Seq125 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Chitinophaga arvensicola
Seq496 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Chitinophaga soli
Seq76 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Chitinophaga NA
Seq410 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Chitinophaga NA
Seq50 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Taibaiella NA
Seq481 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Taibaiella NA
Seq147 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Taibaiella NA
Seq1006 Bacteria Pseudomonadota Gammaproteobacteria Pseudomonadales Moraxellaceae NA NA
Seq727 Bacteria Pseudomonadota Gammaproteobacteria Steroidobacterales Steroidobacteraceae Steroidobacter NA
Seq409 Bacteria Pseudomonadota Gammaproteobacteria Gammaproteobacteria Incertae Sedis Unknown Family Acidibacter NA
Seq471 Bacteria Pseudomonadota Gammaproteobacteria Gammaproteobacteria Incertae Sedis Unknown Family Acidibacter NA
Seq581 Bacteria Pseudomonadota Gammaproteobacteria Gammaproteobacteria Incertae Sedis Unknown Family Acidibacter NA
Seq335 Bacteria Pseudomonadota Gammaproteobacteria Xanthomonadales Rhodanobacteraceae Tahibacter NA
Seq186 Bacteria Pseudomonadota Gammaproteobacteria Xanthomonadales Xanthomonadaceae Xanthomonas NA
Seq606 Bacteria Pseudomonadota Gammaproteobacteria Diplorickettsiales Diplorickettsiaceae NA NA

$darkgreen \(darkgreen\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 29 taxa and 35 samples ] sample_data() Sample Data: [ 35 samples by 4 sample variables ] tax_table() Taxonomy Table: [ 29 taxa by 7 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 29 tips and 28 internal nodes ] refseq() DNAStringSet: [ 29 reference sequences ]

\(darkgreen\)amp ampvis2 object with 5 elements. Summary of OTU table: [1] 35 29 173281 336 15108 4808 4950.89

Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 29(100%) 29(100%) 29(100%) 29(100%) 29(100%) 28(96.55%) 2(6.9%)

Metadata variables: 5 SampleID, Bag, Description, Group, Day

\(darkgreen\)heat \(darkgreen\)heat_rel \(darkgreen\)tree \(darkgreen\)taxa

Kingdom Phylum Class Order Family Genus Species
Seq120 Bacteria Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae Burkholderia telluris
Seq69 Bacteria Pseudomonadota Alphaproteobacteria Azospirillales Inquilinaceae Inquilinus ginsengisoli
Seq4 Bacteria Pseudomonadota Alphaproteobacteria Azospirillales Inquilinaceae Inquilinus NA
Seq370 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Xanthobacteraceae Pseudorhodoplanes NA
Seq9 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Xanthobacteraceae Bradyrhizobium NA
Seq82 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Xanthobacteraceae Bradyrhizobium NA
Seq22 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Xanthobacteraceae Starkeya NA
Seq44 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Rhizobiaceae Allorhizobium NA
Seq10 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Rhizobiaceae Allorhizobium NA
Seq29 Bacteria Myxococcota Polyangia Nannocystales Nannocystaceae Nannocystis NA
Seq364 Bacteria Myxococcota Polyangia Polyangiales Polyangiaceae Labilithrix NA
Seq51 Bacteria Actinobacteriota Actinobacteria Corynebacteriales Mycobacteriaceae Mycobacterium NA
Seq34 Bacteria Actinobacteriota Actinobacteria Streptosporangiales Streptosporangiaceae Herbidospora NA
Seq57 Bacteria Bacillota Bacilli Paenibacillales Paenibacillaceae Paenibacillus NA
Seq25 Bacteria Bacillota Bacilli Bacillales Bacillaceae Terribacillus NA
Seq3 Bacteria Bacillota Bacilli Bacillales Bacillaceae Bacillus NA
Seq5 Bacteria Bacillota Bacilli Bacillales Bacillaceae Bacillus NA
Seq13 Bacteria Bacillota Bacilli Bacillales Planococcaceae NA NA
Seq6 Bacteria Bacillota Bacilli Bacillales Planococcaceae Solibacillus NA
Seq19 Bacteria Bacillota Bacilli Bacillales Planococcaceae Solibacillus NA
Seq17 Bacteria Planctomycetota Planctomycetes Isosphaerales Isosphaeraceae Singulisphaera NA
Seq527 Bacteria Acidobacteriota Acidobacteriae Acidobacteriales Acidobacteriaceae (Subgroup 1) Edaphobacter NA
Seq16 Bacteria Bacteroidota Bacteroidia Cytophagales Microscillaceae Ohtaekwangia NA
Seq7 Bacteria Bacteroidota Bacteroidia Cytophagales Microscillaceae Ohtaekwangia NA
Seq1 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Chitinophaga NA
Seq134 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Chitinophaga NA
Seq81 Bacteria Bacteroidota Bacteroidia Cytophagales Spirosomaceae Dyadobacter NA
Seq15 Bacteria Pseudomonadota Gammaproteobacteria Xanthomonadales Rhodanobacteraceae Luteibacter NA
Seq26 Bacteria Pseudomonadota Gammaproteobacteria Xanthomonadales Xanthomonadaceae Luteimonas NA

$red \(red\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 139 taxa and 35 samples ] sample_data() Sample Data: [ 35 samples by 4 sample variables ] tax_table() Taxonomy Table: [ 139 taxa by 7 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 139 tips and 138 internal nodes ] refseq() DNAStringSet: [ 139 reference sequences ]

\(red\)amp ampvis2 object with 5 elements. Summary of OTU table: [1] 35 139 109849 0 14357 1966 3138.54

Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 139(100%) 139(100%) 139(100%) 139(100%) 139(100%) 86(61.87%) 11(7.91%)

Metadata variables: 5 SampleID, Bag, Description, Group, Day

\(red\)heat \(red\)heat_rel \(red\)tree \(red\)taxa

Kingdom Phylum Class Order Family Genus Species
Seq117 Archaea Thermoproteota Nitrososphaeria Nitrososphaerales Nitrososphaeraceae NA NA
Seq271 Archaea Thermoproteota Nitrososphaeria Nitrososphaerales Nitrososphaeraceae Candidatus Nitrocosmicus NA
Seq61 Archaea Thermoproteota Nitrososphaeria Nitrososphaerales Nitrososphaeraceae NA NA
Seq62 Archaea Thermoproteota Nitrososphaeria Nitrososphaerales Nitrososphaeraceae NA NA
Seq214 Bacteria Pseudomonadota Gammaproteobacteria Diplorickettsiales Diplorickettsiaceae Aquicella NA
Seq758 Bacteria Pseudomonadota Gammaproteobacteria Diplorickettsiales Diplorickettsiaceae Aquicella NA
Seq289 Bacteria Cyanobacteriota Vampirivibrionia Obscuribacterales Obscuribacteraceae NA NA
Seq23 Bacteria Cyanobacteriota Vampirivibrionia Obscuribacterales Obscuribacteraceae NA NA
Seq318 Bacteria Pseudomonadota Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas naasensis
Seq71 Bacteria Pseudomonadota Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas NA
Seq480 Bacteria Pseudomonadota Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas jaspsi
Seq310 Bacteria Pseudomonadota Alphaproteobacteria Dongiales Dongiaceae Dongia NA
Seq262 Bacteria Pseudomonadota Alphaproteobacteria Reyranellales Reyranellaceae Reyranella soli
Seq221 Bacteria Pseudomonadota Alphaproteobacteria Reyranellales Reyranellaceae Reyranella NA
Seq83 Bacteria Pseudomonadota Alphaproteobacteria Reyranellales Reyranellaceae Reyranella NA
Seq321 Bacteria Pseudomonadota Alphaproteobacteria Reyranellales Reyranellaceae Reyranella NA
Seq226 Bacteria Pseudomonadota Alphaproteobacteria Reyranellales Reyranellaceae Reyranella massiliensis
Seq537 Bacteria Pseudomonadota Alphaproteobacteria Rhodospirillales Magnetospiraceae NA NA
Seq712 Bacteria Pseudomonadota Alphaproteobacteria Acetobacterales Acetobacteraceae NA NA
Seq433 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Kaistiaceae Kaistia NA
Seq270 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Xanthobacteraceae Pseudorhodoplanes NA
Seq107 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Xanthobacteraceae Pseudolabrys NA
Seq166 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Xanthobacteraceae Pseudolabrys NA
Seq183 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Xanthobacteraceae Pseudolabrys NA
Seq426 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Xanthobacteraceae NA NA
Seq175 Bacteria Pseudomonadota Alphaproteobacteria Micropepsales Micropepsaceae NA NA
Seq129 Bacteria Pseudomonadota Alphaproteobacteria Micropepsales Micropepsaceae NA NA
Seq157 Bacteria Pseudomonadota Alphaproteobacteria Sphingomonadales Sphingomonadaceae Altererythrobacter NA
Seq137 Bacteria Pseudomonadota Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingobium NA
Seq28 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Devosiaceae Devosia NA
Seq592 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Rhizobiaceae Mesorhizobium NA
Seq711 Bacteria Pseudomonadota Alphaproteobacteria Caulobacterales Caulobacteraceae Phenylobacterium NA
Seq33 Bacteria Pseudomonadota Alphaproteobacteria Caulobacterales Caulobacteraceae Caulobacter NA
Seq225 Bacteria Pseudomonadota Alphaproteobacteria Caulobacterales Hyphomonadaceae Hirschia NA
Seq39 Bacteria Gemmatimonadota Gemmatimonadetes Gemmatimonadales Gemmatimonadaceae NA NA
Seq179 Bacteria Nitrospirota Nitrospiria Nitrospirales Nitrospiraceae Nitrospira japonica
Seq30 Bacteria Actinobacteriota Thermoleophilia Solirubrobacterales Solirubrobacteraceae Conexibacter NA
Seq156 Bacteria Actinobacteriota Thermoleophilia Solirubrobacterales Solirubrobacteraceae Solirubrobacter NA
Seq192 Bacteria Actinobacteriota Thermoleophilia Solirubrobacterales 67-14 NA NA
Seq383 Bacteria Actinobacteriota Thermoleophilia Solirubrobacterales 67-14 NA NA
Seq532 Bacteria Myxococcota Polyangia Haliangiales Haliangiaceae Haliangium NA
Seq504 Bacteria Myxococcota Polyangia Polyangiales BIrii41 NA NA
Seq172 Bacteria Myxococcota Polyangia Polyangiales BIrii41 NA NA
Seq46 Bacteria Myxococcota Polyangia Polyangiales BIrii41 NA NA
Seq493 Bacteria Myxococcota Polyangia Polyangiales BIrii41 NA NA
Seq100 Bacteria Myxococcota Polyangia Polyangiales BIrii41 NA NA
Seq35 Bacteria Myxococcota Polyangia Polyangiales BIrii41 NA NA
Seq399 Bacteria Myxococcota Polyangia Polyangiales Polyangiaceae Minicystis NA
Seq128 Bacteria Myxococcota Polyangia Polyangiales Sandaracinaceae NA NA
Seq274 Bacteria Bdellovibrionota Bdellovibrionia Bdellovibrionales Bdellovibrionaceae Bdellovibrio NA
Seq119 Bacteria Bdellovibrionota Bdellovibrionia Bdellovibrionales Bdellovibrionaceae Bdellovibrio NA
Seq253 Bacteria Dependentiae Babeliae Babeliales UBA12409 NA NA
Seq79 Bacteria Actinobacteriota Actinobacteria Propionibacteriales Propionibacteriaceae Microlunatus NA
Seq416 Bacteria Actinobacteriota Actinobacteria Propionibacteriales Nocardioidaceae Kribbella NA
Seq381 Bacteria Actinobacteriota Acidimicrobiia Microtrichales Ilumatobacteraceae NA NA
Seq173 Bacteria Actinobacteriota Acidimicrobiia Microtrichales Iamiaceae Iamia NA
Seq235 Bacteria Actinobacteriota Acidimicrobiia Microtrichales Iamiaceae Iamia NA
Seq43 Bacteria Actinobacteriota Actinobacteria Micrococcales Microbacteriaceae Galbitalea NA
Seq290 Bacteria Actinobacteriota Actinobacteria Micrococcales Microbacteriaceae Galbitalea NA
Seq693 Bacteria Actinobacteriota Actinobacteria Micromonosporales Micromonosporaceae Luedemannella NA
Seq148 Bacteria Actinobacteriota Actinobacteria Micromonosporales Micromonosporaceae Dactylosporangium NA
Seq357 Bacteria Actinobacteriota Actinobacteria Streptomycetales Streptomycetaceae NA NA
Seq317 Bacteria Bacteroidota Bacteroidia Sphingobacteriales Sphingobacteriaceae Mucilaginibacter NA
Seq230 Bacteria Bacteroidota Bacteroidia Sphingobacteriales Sphingobacteriaceae Mucilaginibacter calamicampi
Seq90 Bacteria Bacteroidota Bacteroidia Sphingobacteriales NS11-12 marine group NA NA
Seq520 Bacteria Bacteroidota Bacteroidia Sphingobacteriales NS11-12 marine group NA NA
Seq145 Bacteria Bacteroidota Bacteroidia Sphingobacteriales env.OPS 17 NA NA
Seq284 Bacteria Bacillota Bacilli Paenibacillales Paenibacillaceae Paenibacillus lautus
Seq258 Bacteria Bacillota Bacilli Paenibacillales Paenibacillaceae Paenibacillus NA
Seq244 Bacteria Bacillota Bacilli Bacillales Planococcaceae Paenisporosarcina NA
Seq85 Bacteria Spirochaetota Spirochaetia Spirochaetales Spirochaetaceae Salinispira NA
Seq185 Bacteria Spirochaetota Spirochaetia Spirochaetales Spirochaetaceae Spirochaeta 2 NA
Seq236 Bacteria Spirochaetota Spirochaetia Spirochaetales Spirochaetaceae Spirochaeta 2 NA
Seq590 Bacteria Patescibacteria Saccharimonadia Saccharimonadales LWQ8 NA NA
Seq238 Bacteria Fibrobacterota Fibrobacteria Fibrobacterales Fibrobacteraceae NA NA
Seq490 Bacteria Chloroflexota Anaerolineae Anaerolineales Anaerolineaceae NA NA
Seq612 Bacteria Chloroflexota Chloroflexia Chloroflexales Roseiflexaceae NA NA
Seq222 Bacteria Chloroflexota Chloroflexia Chloroflexales Roseiflexaceae NA NA
Seq421 Bacteria Chloroflexota Chloroflexia Chloroflexales Roseiflexaceae NA NA
Seq613 Bacteria Chloroflexota Chloroflexia Thermomicrobiales JG30-KF-CM45 NA NA
Seq530 Bacteria Chloroflexota Chloroflexia Thermomicrobiales JG30-KF-CM45 NA NA
Seq388 Bacteria Armatimonadota Fimbriimonadia Fimbriimonadales Fimbriimonadaceae NA NA
Seq576 Bacteria Verrucomicrobiota Verrucomicrobiae Pedosphaerales Pedosphaeraceae NA NA
Seq417 Bacteria Verrucomicrobiota Verrucomicrobiae Opitutales Opitutaceae Lacunisphaera NA
Seq152 Bacteria Verrucomicrobiota Verrucomicrobiae Opitutales Opitutaceae Lacunisphaera limnophila
Seq517 Bacteria Verrucomicrobiota Verrucomicrobiae Opitutales Opitutaceae Opitutus NA
Seq153 Bacteria Verrucomicrobiota Verrucomicrobiae Chthoniobacterales Terrimicrobiaceae Terrimicrobium NA
Seq281 Bacteria Verrucomicrobiota Verrucomicrobiae Chthoniobacterales Terrimicrobiaceae Terrimicrobium NA
Seq528 Bacteria Verrucomicrobiota Verrucomicrobiae Chthoniobacterales Chthoniobacteraceae LD29 NA
Seq659 Bacteria Verrucomicrobiota Verrucomicrobiae Verrucomicrobiales Rubritaleaceae Luteolibacter NA
Seq713 Bacteria Planctomycetota Phycisphaerae Phycisphaerales Phycisphaeraceae NA NA
Seq601 Bacteria Planctomycetota Planctomycetes Pirellulales Pirellulaceae NA NA
Seq839 Bacteria Planctomycetota Planctomycetes Pirellulales Pirellulaceae NA NA
Seq521 Bacteria Planctomycetota Planctomycetes Pirellulales Pirellulaceae NA NA
Seq585 Bacteria Planctomycetota Planctomycetes Pirellulales Pirellulaceae NA NA
Seq195 Bacteria Planctomycetota Planctomycetes Gemmatales Gemmataceae Gemmata NA
Seq404 Bacteria Planctomycetota Planctomycetes Pirellulales Pirellulaceae Pir4 lineage NA
Seq440 Bacteria Planctomycetota Planctomycetes Pirellulales Pirellulaceae Pir4 lineage NA
Seq300 Bacteria Planctomycetota Planctomycetes Planctomycetales Rubinisphaeraceae SH-PL14 NA
Seq292 Bacteria Planctomycetota Planctomycetes Planctomycetales Schlesneriaceae Schlesneria NA
Seq200 Bacteria Acidobacteriota Blastocatellia Blastocatellales Blastocatellaceae NA NA
Seq342 Bacteria Acidobacteriota Vicinamibacteria Vicinamibacterales Vicinamibacteraceae NA NA
Seq434 Bacteria Acidobacteriota Vicinamibacteria Vicinamibacterales Vicinamibacteraceae NA NA
Seq607 Bacteria Bacteroidota Bacteroidia Cytophagales Microscillaceae Ohtaekwangia NA
Seq12 Bacteria Bacteroidota Bacteroidia Cytophagales Microscillaceae Ohtaekwangia NA
Seq196 Bacteria Bacteroidota Bacteroidia Cytophagales Microscillaceae Ohtaekwangia NA
Seq529 Bacteria Bacteroidota Bacteroidia Cytophagales Microscillaceae Ohtaekwangia NA
Seq11 Bacteria Bacteroidota Bacteroidia Cytophagales Microscillaceae NA NA
Seq570 Bacteria Bacteroidota Bacteroidia Sphingobacteriales env.OPS 17 NA NA
Seq178 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Terrimonas NA
Seq92 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Terrimonas NA
Seq326 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Terrimonas NA
Seq768 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Pseudoflavitalea NA
Seq191 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Flavitalea NA
Seq557 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae NA NA
Seq671 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Edaphobaculum NA
Seq459 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Edaphobaculum NA
Seq550 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Edaphobaculum NA
Seq675 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Edaphobaculum NA
Seq547 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Taibaiella NA
Seq834 Bacteria Bacteroidota Bacteroidia Cytophagales Amoebophilaceae Candidatus Amoebophilus NA
Seq760 Bacteria Pseudomonadota Gammaproteobacteria Coxiellales Coxiellaceae Coxiella NA
Seq353 Bacteria Pseudomonadota Gammaproteobacteria Pseudomonadales 211ds20 NA NA
Seq169 Bacteria Pseudomonadota Gammaproteobacteria Steroidobacterales Steroidobacteraceae Steroidobacter flavus
Seq352 Bacteria Pseudomonadota Gammaproteobacteria Steroidobacterales Steroidobacteraceae Steroidobacter NA
Seq506 Bacteria Pseudomonadota Gammaproteobacteria Gammaproteobacteria Incertae Sedis Unknown Family Acidibacter NA
Seq58 Bacteria Pseudomonadota Gammaproteobacteria Gammaproteobacteria Incertae Sedis Unknown Family Acidibacter NA
Seq264 Bacteria Pseudomonadota Gammaproteobacteria Gammaproteobacteria Incertae Sedis Unknown Family Acidibacter NA
Seq231 Bacteria Pseudomonadota Gammaproteobacteria Gammaproteobacteria Incertae Sedis Unknown Family Acidibacter NA
Seq544 Bacteria Pseudomonadota Gammaproteobacteria Xanthomonadales Xanthomonadaceae Luteimonas vadosa
Seq114 Bacteria Pseudomonadota Gammaproteobacteria Xanthomonadales Rhodanobacteraceae Dokdonella ginsengisoli
Seq207 Bacteria Pseudomonadota Gammaproteobacteria Xanthomonadales Rhodanobacteraceae Dokdonella NA
Seq627 Bacteria Pseudomonadota Gammaproteobacteria Salinisphaerales Solimonadaceae NA NA
Seq562 Bacteria Pseudomonadota Gammaproteobacteria Salinisphaerales Solimonadaceae Alkanibacter NA
Seq786 Bacteria Pseudomonadota Gammaproteobacteria Diplorickettsiales Diplorickettsiaceae NA NA
Seq782 Bacteria Pseudomonadota Gammaproteobacteria Diplorickettsiales Diplorickettsiaceae NA NA
Seq339 Bacteria Pseudomonadota Gammaproteobacteria Diplorickettsiales Diplorickettsiaceae NA NA
Seq138 Bacteria Pseudomonadota Gammaproteobacteria Diplorickettsiales Diplorickettsiaceae NA NA
Seq358 Bacteria Pseudomonadota Gammaproteobacteria Legionellales Legionellaceae Legionella NA

$salmon \(salmon\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 71 taxa and 35 samples ] sample_data() Sample Data: [ 35 samples by 4 sample variables ] tax_table() Taxonomy Table: [ 71 taxa by 7 taxonomic ranks ] phy_tree() Phylogenetic Tree: [ 71 tips and 70 internal nodes ] refseq() DNAStringSet: [ 71 reference sequences ]

\(salmon\)amp ampvis2 object with 5 elements. Summary of OTU table: [1] 35 71 115593 29 20719 779 3302.66

Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 71(100%) 71(100%) 71(100%) 71(100%) 71(100%) 67(94.37%) 11(15.49%)

Metadata variables: 5 SampleID, Bag, Description, Group, Day

\(salmon\)heat \(salmon\)heat_rel \(salmon\)tree \(salmon\)taxa

Kingdom Phylum Class Order Family Genus Species
Seq87 Bacteria Pseudomonadota Gammaproteobacteria Burkholderiales Oxalobacteraceae Massilia armeniaca
Seq151 Bacteria Pseudomonadota Gammaproteobacteria Burkholderiales Oxalobacteraceae NA NA
Seq241 Bacteria Pseudomonadota Gammaproteobacteria Burkholderiales Oxalobacteraceae Pseudoduganella NA
Seq170 Bacteria Pseudomonadota Gammaproteobacteria Burkholderiales Oxalobacteraceae Pseudoduganella eburnea
Seq141 Bacteria Pseudomonadota Gammaproteobacteria Burkholderiales Oxalobacteraceae Massilia NA
Seq324 Bacteria Pseudomonadota Gammaproteobacteria Burkholderiales Oxalobacteraceae Massilia NA
Seq208 Bacteria Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae Cupriavidus NA
Seq8 Bacteria Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae Cupriavidus NA
Seq14 Bacteria Pseudomonadota Gammaproteobacteria Burkholderiales Alcaligenaceae Achromobacter NA
Seq55 Bacteria Pseudomonadota Gammaproteobacteria Burkholderiales Comamonadaceae Xylophilus paradoxus
Seq115 Bacteria Pseudomonadota Gammaproteobacteria Burkholderiales Comamonadaceae NA NA
Seq552 Bacteria Pseudomonadota Gammaproteobacteria Burkholderiales Comamonadaceae Rhizobacter NA
Seq412 Bacteria Pseudomonadota Alphaproteobacteria Sphingomonadales Sphingomonadaceae Sphingomonas mucosissima
Seq540 Bacteria Pseudomonadota Alphaproteobacteria Reyranellales Reyranellaceae Reyranella NA
Seq123 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Beijerinckiaceae Microvirga NA
Seq359 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Rhizobiaceae Allorhizobium NA
Seq371 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Devosiaceae Devosia neptuniae
Seq257 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Rhizobiaceae Shinella NA
Seq202 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Rhizobiaceae NA NA
Seq42 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Rhizobiaceae Mycoplana NA
Seq422 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Rhizobiaceae Ensifer NA
Seq21 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Rhizobiaceae Allorhizobium NA
Seq133 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Rhizobiaceae Neorhizobium NA
Seq348 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Rhizobiaceae Allorhizobium azooxidifex
Seq177 Bacteria Pseudomonadota Alphaproteobacteria Rhizobiales Rhizobiaceae NA NA
Seq423 Bacteria Pseudomonadota Alphaproteobacteria Caulobacterales Caulobacteraceae Phenylobacterium mobile
Seq611 Bacteria Myxococcota Polyangia Haliangiales Haliangiaceae Haliangium NA
Seq401 Bacteria Myxococcota Polyangia Polyangiales Polyangiaceae Pajaroellobacter NA
Seq325 Bacteria Actinobacteriota Actinobacteria Micrococcales Promicromonosporaceae Cellulosimicrobium NA
Seq95 Bacteria Actinobacteriota Actinobacteria Micrococcales Microbacteriaceae Microbacterium NA
Seq74 Bacteria Actinobacteriota Actinobacteria Corynebacteriales Mycobacteriaceae Mycobacterium NA
Seq443 Bacteria Actinobacteriota Actinobacteria Glycomycetales Glycomycetaceae Glycomyces NA
Seq78 Bacteria Actinobacteriota Actinobacteria Streptomycetales Streptomycetaceae Streptomyces NA
Seq94 Bacteria Actinobacteriota Actinobacteria Streptomycetales Streptomycetaceae Streptomyces NA
Seq414 Bacteria Actinobacteriota Actinobacteria Streptomycetales Streptomycetaceae Streptomyces NA
Seq181 Bacteria Verrucomicrobiota Verrucomicrobiae Verrucomicrobiales Verrucomicrobiaceae Roseimicrobium NA
Seq218 Bacteria Bacteroidota Bacteroidia Cytophagales Microscillaceae Ohtaekwangia NA
Seq390 Bacteria Bacteroidota Bacteroidia Cytophagales Microscillaceae Ohtaekwangia NA
Seq142 Bacteria Bacteroidota Bacteroidia Cytophagales Microscillaceae Ohtaekwangia NA
Seq255 Bacteria Bacteroidota Bacteroidia Cytophagales Microscillaceae Ohtaekwangia NA
Seq53 Bacteria Bacteroidota Bacteroidia Flavobacteriales Flavobacteriaceae Flavobacterium NA
Seq261 Bacteria Bacteroidota Bacteroidia Flavobacteriales Weeksellaceae Chryseobacterium ginsenosidimutans
Seq233 Bacteria Bacteroidota Bacteroidia Flavobacteriales Weeksellaceae Chryseobacterium NA
Seq446 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Flavitalea NA
Seq97 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Pseudoflavitalea NA
Seq131 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Niastella NA
Seq136 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Niastella NA
Seq167 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Pseudoflavitalea NA
Seq37 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Chitinophaga NA
Seq215 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Chitinophaga NA
Seq248 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Chitinophaga humicola
Seq112 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Chitinophaga NA
Seq89 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Chitinophaga pinensis
Seq2 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Chitinophaga NA
Seq323 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Chitinophaga NA
Seq20 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Chitinophaga NA
Seq445 Bacteria Bacteroidota Bacteroidia Chitinophagales Chitinophagaceae Chitinophaga NA
Seq38 Bacteria Bacteroidota Bacteroidia Cytophagales Spirosomaceae Dyadobacter NA
Seq362 Bacteria Pseudomonadota Gammaproteobacteria Steroidobacterales Steroidobacteraceae Steroidobacter agariperforans
Seq355 Bacteria Pseudomonadota Gammaproteobacteria Xanthomonadales Xanthomonadaceae Pseudoxanthomonas NA
Seq60 Bacteria Pseudomonadota Gammaproteobacteria Xanthomonadales Xanthomonadaceae Stenotrophomonas NA
Seq130 Bacteria Pseudomonadota Gammaproteobacteria Xanthomonadales Xanthomonadaceae Stenotrophomonas NA
Seq305 Bacteria Pseudomonadota Gammaproteobacteria Xanthomonadales Xanthomonadaceae Lysobacter NA
Seq41 Bacteria Pseudomonadota Gammaproteobacteria Xanthomonadales Xanthomonadaceae Lysobacter NA
Seq667 Bacteria Pseudomonadota Gammaproteobacteria Legionellales Legionellaceae Legionella NA
Seq176 Bacteria Pseudomonadota Gammaproteobacteria Pseudomonadales Pseudomonadaceae Pseudomonas NA
Seq24 Bacteria Pseudomonadota Gammaproteobacteria Pseudomonadales Pseudomonadaceae Pseudomonas NA
Seq27 Bacteria Pseudomonadota Gammaproteobacteria Pseudomonadales Pseudomonadaceae Pseudomonas NA
Seq99 Bacteria Pseudomonadota Gammaproteobacteria Pseudomonadales Pseudomonadaceae Pseudomonas NA
Seq56 Bacteria Pseudomonadota Gammaproteobacteria Pseudomonadales Pseudomonadaceae Pseudomonas NA
Seq279 Bacteria Pseudomonadota Gammaproteobacteria Enterobacterales Enterobacteriaceae Klebsiella NA

0.7.3 mdp and other indices for clusters

I can’t do stat mdp support for clusters, but the early clusters(salmon, cyan) have lower mdp than the later ones.

external_empty_dataframe <- data.frame(cluster=factor(),
                                       mdp=numeric(),
                                       richness=numeric(),
                                       stringsAsFactors = FALSE)

for (i in names(l_vst)) {
  ps.in <- l_vst[[i]][["ps"]]
  m <- ps.in@otu_table@.Data %>% 
    colSums() %>% 
    as.data.frame() %>% 
    as.matrix() %>% 
    t()
  
  mdp_index <- picante::mpd(m, cophenetic(ps.in@phy_tree)) 
  ricness_index <- ps.in@tax_table %>% 
    rownames() %>% 
    length()
  d <- data.frame(cluster = i,
                  mdp = mdp_index,
                  richness = ricness_index)
  external_empty_dataframe <- rbind(external_empty_dataframe, d)
}

external_empty_dataframe %>% 
  dplyr::arrange(mdp)
##     cluster      mdp richness
## 1 darkgreen 1.541596       29
## 2    salmon 1.629239       71
## 3      cyan 1.754819       82
## 4       red 1.928198      139
colSums(l_vst$salmon$ps@otu_table@.Data)
##  Seq87 Seq151 Seq241 Seq170 Seq141 Seq324 Seq208   Seq8  Seq14  Seq55 Seq115 
##   1570    975    549    765    703    390    669   9841   6482   2589   1262 
## Seq552 Seq412 Seq540 Seq123 Seq359 Seq371 Seq257 Seq202  Seq42 Seq422  Seq21 
##    195    304    200    669    337    341    516    688   3429    261   5170 
## Seq133 Seq348 Seq177 Seq423 Seq611 Seq401 Seq325  Seq95  Seq74 Seq443  Seq78 
##   1064    363    798    158    160    310    390   1469   1941    249   1822 
##  Seq94 Seq414 Seq181 Seq218 Seq390 Seq142 Seq255  Seq53 Seq261 Seq233 Seq446 
##   1179    301    791    638    317   1038    523   2663    512    572    266 
##  Seq97 Seq131 Seq136 Seq167  Seq37 Seq215 Seq248 Seq112  Seq89   Seq2 Seq323 
##   1457   1047   1054    875   3566    653    532   1311   1524  16855    392 
##  Seq20 Seq445  Seq38 Seq362 Seq355  Seq60 Seq130 Seq305  Seq41 Seq667 Seq176 
##   5184    266   3539    350    357   2372   1105    189   3465    137    799 
##  Seq24  Seq27  Seq99  Seq56 Seq279 
##   4452   4271   1373   2582    457

0.8 Add a metadata

list.files("meta/")
## [1] "cell_realtime_stat.xlsx" "cell_resp_ch_stat.xlsx" 
## [3] "period_legend.xlsx"
realtime.data <- readxl::read_excel("meta/cell_realtime_stat.xlsx")
period.data <- readxl::read_excel("meta/period_legend.xlsx")
resp.data <- readxl::read_excel("meta/cell_resp_ch_stat.xlsx")
realtime.data
## # A tibble: 36 × 29
##    id        day contr…¹ CELL_…² CELL_…³ CELL_…⁴ CELL_…⁵ CELL_…⁶ CELL_…⁷ CELL_…⁸
##    <chr>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
##  1 G01-1-1     3    15.6    33.2    31.1    33.6    36.6    33.8    38.0    33.0
##  2 G01-1-2     3    15.5    33.3    31.3    32.9    36.2    33.8    41      32.7
##  3 G01-1-3     3    15.5    33.6    31.2    32.9    35.8    33.3    38.2    32.9
##  4 G01-2-1     3    17.1    36.0    31.1    33.1    36.9    36.4    38.8    33.8
##  5 G01-2-2     3    17.3    35.6    31.0    32.8    37.8    35.7    40.1    33.6
##  6 G01-2-3     3    17.1    36      31.1    32.8    38.1    35.2    37.6    33.6
##  7 G01-3-1     3    16.4    35.4    31.6    34.9    39.1    35.3    40.9    33.2
##  8 G01-3-2     3    16.5    35.6    31.3    34.4    37.9    34.6    38.9    33.1
##  9 G01-3-3     3    16.4    35.8    31.5    34.4    41.6    35.0    36.5    33.4
## 10 G05-1-1    28    18.2    26.3    29.1    29.4    28.2    28.4    37.4    34.5
## # … with 26 more rows, 19 more variables: CELL_193122 <dbl>, CELL_73229 <dbl>,
## #   CELL_47814 <dbl>, CELL_163125 <dbl>, CELL_73266 <dbl>, CELL_88582 <dbl>,
## #   CELL_63583 <dbl>, CELL_14199 <dbl>, CELL_95616 <dbl>, CELL_63504 <dbl>,
## #   CELL_08643 <chr>, CELL_199599 <dbl>, CELL_01426 <dbl>, CELL_71601 <dbl>,
## #   CELL_45099 <dbl>, CELL_191900 <dbl>, CELL_99463 <dbl>, CELL_74579 <dbl>,
## #   CELL_183403 <dbl>, and abbreviated variable names ¹​contr_16S, ²​CELL_172283,
## #   ³​CELL_203163, ⁴​CELL_83325, ⁵​CELL_188413, ⁶​CELL_109631, ⁷​CELL_188119, …

0.8.1 Realtime data

impute.mean <- function(x) replace(x, is.na(x), mean(x, na.rm = TRUE))

realtime_data <- realtime.data %>% 
  mutate(CELL_08643 = as.numeric(CELL_08643)) %>% 
  group_by(day) %>%
  mutate(nice_cell = impute.mean(CELL_08643)) %>% 
  mutate(day = as.factor(day))
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `CELL_08643 = as.numeric(CELL_08643)`.
## Caused by warning:
## ! NAs introduced by coercion
realtime_zjena <- readxl::read_excel("cellulases_gene_expression(1).xlsx")
geom_mean <-function(x){exp(mean(log(x)))}

realtime_zjena_geom <- realtime_zjena %>% 
  mutate(repeats = paste0(realtime_zjena$week, "-", rep(1:3, 3, each=3))) %>% 
  relocate(repeats, 1) %>% 
  group_by(repeats) %>% 
  summarise_if(is.numeric, geom_mean) %>% 
  mutate(week = as.factor(week)) %>% 
  arrange(week)

UPGMA cellulase(distance - bray).
I don’t really know what that means. Depending on the distance, you get very different clusters.

realtime_matrix <- realtime_zjena_geom %>% 
  column_to_rownames("repeats") %>% 
  select_if(is.numeric) %>% 
  as.matrix()

hcl <- hclust(vegan::vegdist(t(realtime_matrix), method="bray"), "average")
plot(hcl)

And here is the Euclidean distance.
Nothing to compare with the previous picture

hcl <- hclust(vegan::vegdist(t(realtime_matrix), method="euclidian"), "average")
plot(hcl)

Cellulases corplot.
Correlations only

m = cor(realtime_matrix)
corrplot::corrplot(m)

Only significant correlations . In general, there are clusters, but the correlations are not significant (with some exceptions)

cor_test_mat <- psych::corr.test(realtime_matrix)$p
corrplot::corrplot(m, p.mat = cor_test_mat, method = 'circle', type = 'lower', insig='blank',
         order = 'AOE', diag = FALSE)$corrPos -> p1

# text(p1$x, p1$y, round(p1$corr, 2))

The same, but the matrix is already logarithmased

cor_test_mat <- psych::corr.test(log(realtime_matrix))
corrplot::corrplot(cor_test_mat$r, p.mat = cor_test_mat$p, method = 'circle', type = 'lower', insig='blank',
         order = 'AOE', diag = FALSE)$corrPos -> p1

# text(p1$x, p1$y, round(p1$corr, 2))

add phylogenic tree (mafft - iqtree(ML))

realtime_tree <- ape::read.tree("al_chz.fasta.contree")
plot(realtime_tree)

library(tidyverse)

realtime_data %>% 
  select(c("day", "id", "contr_16S", "CELL_172283")) %>% 
  mutate(bio_repl = gsub("-[1-3]$", "", id)) %>% 
  group_by(bio_repl, day) %>% 
  summarise(contr_tmean = mean(contr_16S),
            data_tmean = mean(CELL_172283)) %>% 
  mutate(dCt = data_tmean - contr_tmean)
## `summarise()` has grouped output by 'bio_repl'. You can override using the
## `.groups` argument.
## # A tibble: 12 × 5
## # Groups:   bio_repl [12]
##    bio_repl day   contr_tmean data_tmean   dCt
##    <chr>    <fct>       <dbl>      <dbl> <dbl>
##  1 G01-1    3            15.5       33.4 17.8 
##  2 G01-2    3            17.2       35.9 18.7 
##  3 G01-3    3            16.4       35.6 19.2 
##  4 G05-1    28           18.4       26.2  7.81
##  5 G05-2    28           16.5       26.2  9.73
##  6 G05-3    28           18.5       32.0 13.6 
##  7 G10-1    91           17.4       26.2  8.72
##  8 G10-2    91           16.3       24.1  7.79
##  9 G10-3    91           17.5       26.3  8.81
## 10 G14-1    161          18.0       28.1 10.1 
## 11 G14-2    161          16.5       27.0 10.6 
## 12 G14-3    161          16.9       25.4  8.41
resp_data <- resp.data %>% 
  group_by(day) %>%
  mutate(control = impute.mean(control)) %>% 
  mutate(straw = impute.mean(straw)) 

resp_data
## # A tibble: 164 × 3
## # Groups:   day [27]
##      day control straw
##    <dbl>   <dbl> <dbl>
##  1     0    250   250 
##  2     3    285.  723.
##  3     3    263.  832.
##  4     3    131.  657.
##  5     3    307.  525.
##  6     3    241.  744.
##  7     3    245.  701.
##  8     3    245.  657.
##  9     7    274.  690.
## 10     7    296.  690.
## # … with 154 more rows

0.8.2 Knee plot

Soil resperation - median(straw)/median(control)

You can use a third-party package (KneeArrower) to understand where the knee_plot breaks. For the elimination of repetitions let’s take the median.   I don’t know how this package works (looks for a derivative, but I don’t know how it smoothes it out, it’s math), but it says that the breaking is more likely at 60-80 days.

(I corrected the errors - got closer to your data)

https://github.com/agentlans/KneeArrower - here’s where you can read about math

# resp_data %>% 
#   filter(!day == 0) %>% 
#   group_by(day) %>% 
#   summarise(median_control = median(control),
#             median_straw = median(straw)) %>% 
#   mutate(rel = median_straw/median_control) %>% 
#   ggplot() +
#   geom_point(aes(x = day, y = rel))

resp_median <- resp_data %>% 
  filter(!day == 0) %>% 
  group_by(day) %>% 
  summarise(median_control = median(control),
            median_straw = median(straw)) %>% 
  mutate(rel = median_straw/median_control)


thresholds <- c(0.25, 0.5, 0.75, 1)

# Find cutoff points at each threshold
cutoff.points <- lapply(thresholds, function(i) {
  KneeArrower::findCutoff(resp_median$day, resp_median$rel, method="first", i)
})

x.coord <- sapply(cutoff.points, function(p) p$x)
y.coord <- sapply(cutoff.points, function(p) p$y)

# Plot the cutoff points on the scatterplot
plot(resp_median$day, resp_median$rel, pch=20, col="gray")
points(x.coord, y.coord, col="red", pch=20)
text(x.coord, y.coord, labels=thresholds, pos=4, col="red")

period_data <- period.data %>%
  mutate(bag_id = as.factor(bag_id))

period_data
## # A tibble: 15 × 2
##    bag_id   day
##    <fct>  <dbl>
##  1 1          3
##  2 2          7
##  3 3         14
##  4 4         21
##  5 5         28
##  6 6         35
##  7 7         49
##  8 8         63
##  9 9         77
## 10 10        91
## 11 11       105
## 12 12       119
## 13 13       140
## 14 14       161
## 15 15       182

0.8.3 alpha/respiration.

Well, I don’t really understand what to do next - attach the clusters to this picture?

resp_median_bags <- resp_median %>% 
  left_join(period.data, by="day") %>% 
  mutate(bag_id = as.factor(bag_id))

ps.f.r <- rarefy_even_depth(ps.f)

estimate_richness(ps.f.r) %>% 
  rownames_to_column("ID") %>% 
  mutate(bag_id = as.factor(
    as.numeric(
      gsub("\\..+","",
           gsub("straw\\.16s\\.D","", ID)
           )
      )
    )
    ) %>% 
  group_by(bag_id) %>% 
  summarise(Observed = mean(Observed),
            Shannon = mean(Shannon),
            InvSimpson = mean(InvSimpson)) %>% 
  left_join(period_data, by="bag_id") %>% 
  left_join(resp_median_bags, by="bag_id") %>% 
  mutate(
    Observed = scale(Observed),
    Shannon = scale(Shannon),
    InvSimpson = scale(InvSimpson),
    Respiration = scale(rel)
  ) %>% 
  select(c(bag_id, Observed, Shannon, InvSimpson, Respiration)) %>% 
  pivot_longer(c("Observed", "Shannon", "InvSimpson", "Respiration")) %>% 
  mutate(metric = factor(name, levels = c("Observed", "Shannon", "InvSimpson", "Respiration"))) %>% 
  ggplot(aes(y = value, x = bag_id, group = name)) +
  geom_line(aes(color = name),
            alpha = 0.8) +
  theme_bw() +
    scale_x_discrete(labels=(c("1"="3",
                                "3"="14",
                                "5"="28",
                                "7"="49" ,
                                "8"="63",
                                "10"="91",
                                "12"="119",
                                "13"="140",
                                "14"="161",
                                "15"="182")
         )) +
  labs(
    x = "Days",
    y = "z-scaled mean values",
    colour = ""
   )

Correlation is negative, weak, not particularly reliable, and only Pearson (which works for a normal distribution (we have a time series, we need a spearman for good measure))

alpha_resp <- phyloseq::estimate_richness(ps.f.r) %>% 
  rownames_to_column("ID") %>% 
  mutate(bag_id = as.factor(
    as.numeric(
      gsub("\\..+","",
           gsub("straw\\.16s\\.D","", ID)
           )
      )
    )
    ) %>% 
  group_by(bag_id) %>% 
  summarise(Observed = mean(Observed),
            Shannon = mean(Shannon),
            InvSimpson = mean(InvSimpson)) %>% 
  left_join(period_data, by="bag_id") %>% 
  left_join(resp_median_bags, by="bag_id") %>% 
  mutate(
    Observed_scaled = scale(Observed),
    Shannon_scaled = scale(Shannon),
    InvSimpson_scaled = scale(InvSimpson),
    Respiration_scaled = scale(rel)
  ) %>% 
  select(c(bag_id, Observed_scaled, Shannon_scaled, InvSimpson_scaled, Respiration_scaled))

cor.test(alpha_resp$Observed_scaled, alpha_resp$Respiration_scaled, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  alpha_resp$Observed_scaled and alpha_resp$Respiration_scaled
## S = 240, p-value = 0.1909
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.4545455
cor.test(alpha_resp$Shannon_scaled, alpha_resp$Respiration_scaled, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  alpha_resp$Shannon_scaled and alpha_resp$Respiration_scaled
## S = 242, p-value = 0.1782
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.4666667
cor.test(alpha_resp$InvSimpson_scaled, alpha_resp$Respiration_scaled, method = "spearman")
## 
##  Spearman's rank correlation rho
## 
## data:  alpha_resp$InvSimpson_scaled and alpha_resp$Respiration_scaled
## S = 242, p-value = 0.1782
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.4666667
cor.test(alpha_resp$Observed_scaled, alpha_resp$Respiration_scaled, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  alpha_resp$Observed_scaled and alpha_resp$Respiration_scaled
## t = -2.5961, df = 8, p-value = 0.03181
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.91588600 -0.08110917
## sample estimates:
##        cor 
## -0.6762028
cor.test(alpha_resp$Shannon_scaled, alpha_resp$Respiration_scaled, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  alpha_resp$Shannon_scaled and alpha_resp$Respiration_scaled
## t = -2.858, df = 8, p-value = 0.02122
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9259959 -0.1468753
## sample estimates:
##        cor 
## -0.7107719
cor.test(alpha_resp$InvSimpson_scaled, alpha_resp$Respiration_scaled, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  alpha_resp$InvSimpson_scaled and alpha_resp$Respiration_scaled
## t = -2.2456, df = 8, p-value = 0.05494
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.89933340  0.01287072
## sample estimates:
##        cor 
## -0.6217944

0.9 Export all phyloseq object in one table

In the table there is a column cluster - “exl” in it means that this is the part of the dataset that did not go to WGCNA - majors in single samples(OMs).

ps.m <- phyloseq::psmelt(ps.f)

ps.m <- ps.m %>% 
  mutate_if(is.character, as.factor) 

ps.data.out <- ps.m %>%
  select(-Group) %>% 
  pivot_wider(names_from = c(Day, Description, Sample), values_from = Abundance, values_fill = 0) 

#create empty dataframe with columnnames
external_empty_dataframe <- data.frame(OTU=factor(), cluster=factor(), stringsAsFactors = FALSE)

for (i in names(l_vst)) {
  a <- taxa_names(l_vst[[i]][["ps"]])
  b <- rep(i, length(a))
  d <- data.frame(OTU = as.factor(a),
                  cluster = as.factor(b))
  external_empty_dataframe <- rbind(external_empty_dataframe, d)
}

clusters.otu.df <- external_empty_dataframe
# add exl taxa -- taxa "exl"
d <- data.frame(OTU = as.factor(ps.exl.taxa),
                cluster = as.factor(rep("exl", length(ps.exl.taxa))))

clusters.otu.df <- rbind(clusters.otu.df, d)
ps.data.out.exl <- left_join(clusters.otu.df, ps.data.out, by="OTU")
## Warning in left_join(clusters.otu.df, ps.data.out, by = "OTU"): Each row in `x` is expected to match at most 1 row in `y`.
## ℹ Row 1 of `x` matches multiple rows.
## ℹ If multiple matches are expected, set `multiple = "all"` to silence this
##   warning.
write.table(ps.data.out.exl, file = "ps.data.out.tsv", sep = "\t")
ps.data.out.exl

0.9.1 Several pictures on the dynamics of representation of individual phyla

Bdellovibrionota - predators, an indicator of a developed community

ps.m %>% 
  filter(Phylum == "Bdellovibrionota") %>% 
  group_by(Description, Day) %>% 
  summarise(Bs = sum(Abundance)) %>% 
  ggplot() +
  geom_boxplot(aes(x = Day, y = Bs)) +
  theme_bw()
## `summarise()` has grouped output by 'Description'. You can override using the
## `.groups` argument.

Myxococcota - sort of the same, but as we know can be cellulotic(facultative predators)

ps.m %>% 
  filter(Phylum == "Myxococcota") %>% 
  group_by(Description, Day) %>% 
  summarise(Bs = sum(Abundance)) %>% 
  ggplot() +
  geom_boxplot(aes(x = Day, y = Bs)) +
  theme_bw()
## `summarise()` has grouped output by 'Description'. You can override using the
## `.groups` argument.

Археи появляются тоже на поздних стадиях - вообще я бы хотел бы опять развить тему важности азотного метаболизма на поздних стадиях разложения.
Оч хочется метагеном, но не этот.
Кроме того хочется отметить, что эти минорные группы возникают на D12

ps.m %>% 
  filter(Phylum == "Crenarchaeota") %>% 
  group_by(Description, Day) %>% 
  summarise(Bs = sum(Abundance)) %>% 
  ggplot() +
  geom_boxplot(aes(x = Day, y = Bs)) +
  theme_bw()
## `summarise()` has grouped output by 'Description'. You can override using the
## `.groups` argument.

Gammaproteobacteria - they are a third of the cluster blue

ps.m %>% 
  filter(Class == "Gammaproteobacteria") %>% 
  group_by(Description, Day) %>% 
  summarise(Bs = sum(Abundance)) %>% 
  ggplot() +
  geom_boxplot(aes(x = Day, y = Bs)) +
  theme_bw()
## `summarise()` has grouped output by 'Description'. You can override using the
## `.groups` argument.

0.9.2 Boxplots by phyla

All phyla - the representation is logarithmic on the basis of 2
Separate points are sums of absolute values of reads by days
The sums are logarithmic, but not the individual phylotypes

#select only major phylums
top_phylum <- ps.m %>% 
  count(Phylum) %>% 
  arrange(desc(n)) %>% 
  top_n(10) %>% 
  pull(Phylum)
## Selecting by n
ps.m %>% 
  filter(Phylum %in% top_phylum) %>% 
  mutate(
    Phylum = as.character(Phylum),
    Class = as.character(Class),
    phylum = ifelse(Phylum == "Proteobacteria", Class, Phylum)
  ) %>% 
  group_by(Description, Day, phylum) %>% 
  filter(!is.na(phylum)) %>% 
  summarise(Bs = log2(sum(Abundance))) %>% 
  ggplot(aes(x = Day, y = Bs)) +
  geom_boxplot(fill="#4DBBD5B2", alpha=0.4) +
  theme_bw() +
  facet_wrap(~ phylum)
## `summarise()` has grouped output by 'Description', 'Day'. You can override
## using the `.groups` argument.
## Warning: Removed 54 rows containing non-finite values (`stat_boxplot()`).

0.9.3 resperation and phylotypes correlations

little_res <- resp_median_bags %>% 
  filter(day %in% ps.varstab@sam_data$Day)

ps.varstab.merged <- phyloseq::merge_samples(ps.varstab, "Day")
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'tidytree'
otus_cor <- ps.varstab.merged@otu_table %>% 
  as.data.frame()

s1 <- sapply(colnames(otus_cor), function(x) {
  cor <- cor.test(otus_cor[[x]], little_res$rel, method="spearman", exact = FALSE)
  c1 <- cor$p.value
  c2 <- cor$estimate %>% unname()
  return(list("p.value"=c1, "rho"=c2))})

s1 <- s1 %>% 
  t() %>% 
  as.data.frame() %>% 
  unnest(cols = c(p.value, rho), .id = "id") %>% 
  column_to_rownames("id") %>% 
  mutate(p.adj = p.adjust(p.value, method = "BH", n = length(p.value)))


taxa <- as.data.frame(ps.varstab@tax_table@.Data) %>% 
  rownames_to_column("ID")

s1 %>%
  rownames_to_column("ID") %>% 
  left_join(taxa) %>% 
  arrange(rho) %>%
  filter(Phylum %in% c("Proteobacteria", "Acidobacteriota", "Actinobacteriota","Chloroflexi", "Planctomycetota", "Bacteroidota", "Verrucomicrobiota", "Firmicutes")) %>% 
  mutate(Rank = ifelse(Phylum == "Proteobacteria", Class, Phylum)) %>%
  ggplot(aes(x=rho, color=Rank, fill=Rank)) +
  geom_histogram(alpha=0.3, bins = 30) +
  theme_minimal() +
  facet_wrap(~ Rank) +
  theme(legend.position="none")
## Joining with `by = join_by(ID)`

s1.otus <- s1 %>% 
  rownames_to_column("OTU")
ps.data.out.exl %>% 
  select_if(Negate(is.integer)) %>% 
  select(-Bag) %>% 
  distinct() %>% 
  left_join(s1.otus, by="OTU") %>%
  filter(rho >= 0) %>%
  arrange(p.adj) %>% 
  head()
##     OTU cluster  Kingdom         Phylum               Class           Order
## 1 Seq24  salmon Bacteria Pseudomonadota Gammaproteobacteria Pseudomonadales
## 2 Seq24  salmon Bacteria Pseudomonadota Gammaproteobacteria Pseudomonadales
## 3 Seq24  salmon Bacteria Pseudomonadota Gammaproteobacteria Pseudomonadales
## 4 Seq24  salmon Bacteria Pseudomonadota Gammaproteobacteria Pseudomonadales
## 5 Seq24  salmon Bacteria Pseudomonadota Gammaproteobacteria Pseudomonadales
## 6 Seq24  salmon Bacteria Pseudomonadota Gammaproteobacteria Pseudomonadales
##             Family       Genus Species 91_D10_5_straw.16s.D10.5
## 1 Pseudomonadaceae Pseudomonas    <NA>                        0
## 2 Pseudomonadaceae Pseudomonas    <NA>                        0
## 3 Pseudomonadaceae Pseudomonas    <NA>                        0
## 4 Pseudomonadaceae Pseudomonas    <NA>                        0
## 5 Pseudomonadaceae Pseudomonas    <NA>                        0
## 6 Pseudomonadaceae Pseudomonas    <NA>                        0
##   14_D03_3_straw.16s.D03.3 49_D07_3_straw.16s.D07.3 28_D05_5_straw.16s.D05.5
## 1                     1134                        0                        0
## 2                        0                        0                        0
## 3                        0                        0                      337
## 4                        0                        0                        0
## 5                        0                        0                        0
## 6                        0                        0                        0
##   14_D03_2_straw.16s.D03.2 28_D05_4_straw.16s.D05.4 63_D08_3_straw.16s.D08.3
## 1                     1150                        0                        0
## 2                        0                        0                        0
## 3                        0                        0                        0
## 4                        0                        0                        0
## 5                        0                        0                        0
## 6                        0                        0                        0
##   28_D05_2_straw.16s.D05.2 63_D08_1_straw.16s.D08.1 91_D10_1_straw.16s.D10.1
## 1                        0                        0                        0
## 2                        0                        0                        0
## 3                      367                        0                        0
## 4                        0                        0                        0
## 5                        0                        0                        0
## 6                        0                        0                        0
##   3_D01_2_straw.16s.D01.2 63_D08_2_straw.16s.D08.2 140_D13_2_straw.16s.D13.2
## 1                       0                        0                         0
## 2                      71                        0                         0
## 3                       0                        0                         0
## 4                       0                        0                         0
## 5                       0                        0                         0
## 6                       0                        0                         0
##   119_D12_1_straw.16s.D12.1 91_D10_4_straw.16s.D10.4 49_D07_2_straw.16s.D07.2
## 1                         0                        0                        0
## 2                         0                        0                        0
## 3                         0                        0                        0
## 4                         0                        0                        0
## 5                         0                        0                       53
## 6                         0                        0                        0
##   119_D12_3_straw.16s.D12.3 182_D15_2_straw.16s.D15.2 161_D14_2_straw.16s.D14.2
## 1                         0                         0                         0
## 2                         0                         0                         0
## 3                         0                         0                         0
## 4                         0                         0                         0
## 5                         0                         0                         0
## 6                         0                         0                         0
##   119_D12_2_straw.16s.D12.2 3_D01_3_straw.16s.D01.3 28_D05_1_straw.16s.D05.1
## 1                         0                       0                        0
## 2                         0                     394                        0
## 3                         0                       0                        0
## 4                         0                       0                        0
## 5                         0                       0                        0
## 6                         0                       0                        0
##   91_D10_2_straw.16s.D10.2 3_D01_1_straw.16s.D01.1 91_D10_3_straw.16s.D10.3
## 1                        0                       0                        0
## 2                        0                     370                        0
## 3                        0                       0                        0
## 4                        0                       0                        0
## 5                        0                       0                        0
## 6                        0                       0                        0
##   161_D14_3_straw.16s.D14.3 28_D05_3_straw.16s.D05.3 161_D14_5_straw.16s.D14.5
## 1                         0                        0                         0
## 2                         0                        0                         0
## 3                         0                      377                         0
## 4                         0                        0                         0
## 5                         0                        0                         0
## 6                         0                        0                         0
##   182_D15_3_straw.16s.D15.3 140_D13_1_straw.16s.D13.1 140_D13_3_straw.16s.D13.3
## 1                         0                         0                         0
## 2                         0                         0                         0
## 3                         0                         0                         0
## 4                         0                         0                         0
## 5                         0                         0                         0
## 6                         0                         0                         0
##   161_D14_1_straw.16s.D14.1 14_D03_1_straw.16s.D03.1 161_D14_4_straw.16s.D14.4
## 1                         0                      117                         0
## 2                         0                        0                         0
## 3                         0                        0                         0
## 4                         0                        0                         0
## 5                         0                        0                         0
## 6                         0                        0                         0
##   182_D15_1_straw.16s.D15.1     p.value       rho      p.adj
## 1                         0 0.005446342 0.8000947 0.09957456
## 2                         0 0.005446342 0.8000947 0.09957456
## 3                         0 0.005446342 0.8000947 0.09957456
## 4                        82 0.005446342 0.8000947 0.09957456
## 5                         0 0.005446342 0.8000947 0.09957456
## 6                         0 0.005446342 0.8000947 0.09957456

0.10 Additions

Part added after rounds of peer-review. Thank you very match the reviewer two, you are great!

realtime.stat <- readxl::read_excel("cell_realtime_stat_ph.xlsx")
realtime.names <- readxl::read_excel("cell_realtime_stat_gh_names.xlsx")

realtime.names
## # A tibble: 23 × 4
##    Cellulase   Taxa_Genus   GH_Family Taxa_Family     
##    <chr>       <chr>        <chr>     <chr>           
##  1 CELL_172842 Devosia      GH39      Devosiaceae     
##  2 CELL_203163 Chitinophaga GH43      Chitinophagaceae
##  3 CELL_83325  Chitinophaga GH43      Chitinophagaceae
##  4 CELL_188413 Chitinophaga GH31      Chitinophagaceae
##  5 CELL_109631 Chitinophaga GH5       Chitinophagaceae
##  6 CELL_193122 Chitinophaga GH3       Chitinophagaceae
##  7 CELL_73229  Chitinophaga GH31      Chitinophagaceae
##  8 CELL_47814  Ohtaekwangia GH3       Fulvivirgaceae  
##  9 CELL_163127 Ohtaekwangia GH5       Fulvivirgaceae  
## 10 CELL_73266  Chitinophaga GH28      Chitinophagaceae
## # … with 13 more rows
library(ggpubr)
library(rstatix)
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
## 
##     filter
real_long <- pivot_longer(realtime.stat, colnames(realtime.stat)[2:length(colnames(realtime.stat))]) %>%
  mutate(phase = factor(phase, levels = c("early", "middle", "late")))

stat.test <- real_long %>%
  group_by(name) %>%
  tukey_hsd(value ~ phase) %>%
  add_xy_position(x = "name",
                  dodge = 0.8,
                  step.increase = 0.05)

real_long %>% 
  ggplot(aes(x=name, y=value, col=phase)) +
    geom_boxplot() +
    geom_point(size=0.4, alpha=0.9, position=position_dodge(width=0.75),aes(group=phase)) +
    theme_light() + 
    scale_color_brewer(palette="Dark2") +
    theme(axis.text.x = element_text(angle = 90),
          axis.title.x=element_blank()) + 
  labs( title = "Realtime of selected GH") +
  ylab("log2foldchange") +
  stat_pvalue_manual(
    stat.test,
    label = "p.adj.signif", 
    tip.length = 0.005,
    hide.ns = TRUE,
    bracket.nudge.y = -0.5
  )

realtime.big <- readxl::read_excel("straw_rt_stat_gr.xlsx")
## New names:
## • `` -> `...1`
realtime.big <- realtime.big %>% 
  group_by(day) %>% 
  mutate_at(vars(FUN), 
            ~replace_na(., 
                        mean(., na.rm = TRUE)))

realtime.big.long <- realtime.big %>% 
  filter(phase != "no") %>% 
  pivot_longer(colnames(realtime.big)[4:length(colnames(realtime.big))]) %>% 
  mutate(phase = factor(phase, levels = c("early", "middle", "late"))) %>% 
  mutate(value = log2(value)) 

stat.test <- realtime.big.long %>%
  group_by(name) %>%
  tukey_hsd(value ~ phase) %>%
  add_xy_position(x = "name", dodge = 0.5)

realtime.big.long %>% 
  ggplot(aes(x=name, y=value, col=phase)) +
    geom_boxplot() +
    geom_point(size=0.4, alpha=0.9, position=position_dodge(width=0.75),aes(group=phase)) +
    theme_light() + 
    scale_color_brewer(palette="Dark2") +
    theme(axis.text.x = element_text(angle = 90),
          axis.title.x=element_blank()) + 
  labs(
  title = "Realtime"
) +
ylab("operons per gram, log2") +
  stat_pvalue_manual(
  stat.test,  label = "{p.adj} {p.adj.signif}", 
  tip.length = 0.005, hide.ns = TRUE
  )

gh_taxa <- readr::read_tsv("kraken_gh_full_contig_207_corrected.csv") %>% 
  mutate_if(., is.character, as.factor)
## Rows: 23 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (8): Cellulase, Taxa_Genus, GH_Family, Taxa_Family, Contig_name, Family_...
## dbl (1): Contig_identifier
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
gh_taxa %>% 
  DT::datatable(caption = "Kraken2 based taxonomy for GH containing contigs")

Permanova based on GH specific RealTime data

real_matrix <- realtime.stat %>% 
  select(-phase) %>% 
  as.matrix() %>% 
  t()

gh_meta <- gh_taxa %>% 
  column_to_rownames("Cellulase") %>% 
  select(GH_Family, Genus_new, Family_new) 

meta <- gh_taxa %>% 
  select(GH_Family, Genus_new, Family_new) %>% 
  t() %>% 
  as.data.frame() %>% 
  setNames(gh_taxa$Cellulase)
  
d <- dist(real_matrix, method = "euclidian")

table1 <- phyloseq::distance(ps.f, "bray")
table2 <- as(sample_data(ps.f@sam_data), "data.frame")

list(
  vegan::adonis2(d ~ GH_Family + Genus_new + Family_new, data = gh_meta),
  vegan::adonis2(d ~ Genus_new  + GH_Family, data = gh_meta),
  vegan::adonis2(d ~ GH_Family + Genus_new, data = gh_meta),
  vegan::adonis2(d ~ GH_Family*Genus_new, data = gh_meta)
)
## [[1]]
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = d ~ GH_Family + Genus_new + Family_new, data = gh_meta)
##            Df SumOfSqs      R2      F Pr(>F)   
## GH_Family   7   3479.3 0.41290 3.1769  0.013 * 
## Genus_new   6   2712.2 0.32186 2.8892  0.018 * 
## Family_new  1    983.4 0.11670 6.2854  0.007 **
## Residual    8   1251.7 0.14854                 
## Total      22   8426.6 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[2]]
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = d ~ Genus_new + GH_Family, data = gh_meta)
##           Df SumOfSqs      R2      F Pr(>F)  
## Genus_new  6   4625.9 0.54896 3.1045  0.014 *
## GH_Family  7   1565.6 0.18580 0.9006  0.583  
## Residual   9   2235.1 0.26524                
## Total     22   8426.6 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[3]]
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = d ~ GH_Family + Genus_new, data = gh_meta)
##           Df SumOfSqs      R2      F Pr(>F)  
## GH_Family  7   3479.3 0.41290 2.0015  0.084 .
## Genus_new  6   2712.2 0.32186 1.8202  0.103  
## Residual   9   2235.1 0.26524                
## Total     22   8426.6 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [[4]]
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = d ~ GH_Family * Genus_new, data = gh_meta)
##                     Df SumOfSqs      R2      F Pr(>F)  
## GH_Family            7   3479.3 0.41290 2.2448  0.070 .
## Genus_new            6   2712.2 0.32186 2.0415  0.080 .
## GH_Family:Genus_new  3    906.6 0.10758 1.3648  0.292  
## Residual             6   1328.5 0.15766                
## Total               22   8426.6 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

0.11 ITS

Same analysis, but much quicker and dirty.

library(DESeq2)

ps_its <- readRDS("../../its/d2/ps_its")

sample.data <- ps_its@sam_data %>% 
  data.frame() %>% 
  mutate(Group = if_else(Day %in% c("D01", "D03", "D05"), "early", 
                         if_else(Day %in% c("D07", "D08","D10"), "middle", "late"))) %>% 
  mutate(Group = factor(Group, levels=c("early", "middle","late"))) %>% 
  mutate(Day = Day %>% 
           forcats::fct_recode( "3" = "D01",
                                "14" = "D03",
                                "28" = "D05",
                                "49" = "D07" ,
                                "63" = "D08",
                                "91" = "D10",
                                "119" = "D12",
                                "140" = "D13",
                                "161" = "D14",
                                "182" = "D15")
         ) %>% 
  phyloseq::sample_data()

sample_data(ps_its) <- sample.data

amp_its <- phyloseq_to_ampvis2(ps_its)

ps.its.inall <- phyloseq::filter_taxa(ps_its, function(x) sum(x > 10) > (0.1*length(x)), TRUE)
ps.its.inall <- prune_taxa(taxa_sums(ps.its.inall) > 0, ps.its.inall)
Biostrings::writeXStringSet(ps.its.inall@refseq, file="ref_inall.fasta")
Biostrings::writeXStringSet(ps_its@refseq, file="ref.fasta")

tree <- ape::read.tree("../../its/d2/al.fasta.contree")

ps.its.inall@phy_tree <- tree
ps.its.inall  <- prune_taxa(taxa_sums(ps.its.inall) > 0, ps.its.inall)


ps.its.exl  <- phyloseq::filter_taxa(ps_its, function(x) sum(x > 10) < (0.1*length(x)), TRUE)
ps.its.exl  <- prune_taxa(taxa_sums(ps.its.exl) > 100, ps.its.exl)
ps.its.exl.taxa <- taxa_names(ps.its.exl)

amp_its_inall <- phyloseq_to_ampvis2(ps.its.inall)

Base statistics for its

amp_its
## ampvis2 object with 4 elements. 
## Summary of OTU table:
## [1]       36     1295   815134    13533    34994    22657 22642.61
## 
## Assigned taxonomy:
##     Kingdom      Phylum       Class       Order      Family       Genus 
##    554(43%) 453(34.98%) 362(27.95%) 335(25.87%) 321(24.79%) 306(23.63%) 
##     Species 
## 207(15.98%) 
## 
## Metadata variables: 4 
##  SampleID, Day, Description, Group

Major ITS phylotypes

amp_heatmap(amp_its,
            tax_show = 40,
            group_by = "Day",
            tax_aggregate = "OTU",
            tax_add = "Genus",
            normalise=TRUE,
            showRemainingTaxa = TRUE)
## Warning: Transformation introduced infinite values in discrete y-axis

All by the phylums

amp_heatmap(amp_its,
            tax_show = 7,
            group_by = "Day",
            tax_aggregate = "Phylum",
            tax_add = "Kingdom",
            normalise=TRUE,
            showRemainingTaxa = TRUE)
## Warning: Transformation introduced infinite values in discrete y-axis

There are going to WGCNA pipeline(116 phylotypes)

amp_heatmap(amp_its_inall,
            tax_show = 40 ,
            group_by = "Day",
            tax_aggregate = "OTU",
            tax_add = "Genus",
            normalise=TRUE,
            showRemainingTaxa = TRUE)
## Warning: Transformation introduced infinite values in discrete y-axis

beta_custom_norm_NMDS_elli_w(ps_its, Group="Day", Color="Day")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1834317 
## Run 1 stress 0.1834317 
## ... New best solution
## ... Procrustes: rmse 0.0000328311  max resid 0.00009161747 
## ... Similar to previous best
## Run 2 stress 0.2097497 
## Run 3 stress 0.1841584 
## Run 4 stress 0.1841585 
## Run 5 stress 0.1834317 
## ... New best solution
## ... Procrustes: rmse 0.00001495322  max resid 0.0000474614 
## ... Similar to previous best
## Run 6 stress 0.2231478 
## Run 7 stress 0.1834317 
## ... Procrustes: rmse 0.00002320338  max resid 0.00005725448 
## ... Similar to previous best
## Run 8 stress 0.1841586 
## Run 9 stress 0.1841584 
## Run 10 stress 0.1834317 
## ... Procrustes: rmse 0.00003210621  max resid 0.00008949888 
## ... Similar to previous best
## Run 11 stress 0.1841584 
## Run 12 stress 0.1834317 
## ... Procrustes: rmse 0.000003284702  max resid 0.00001199426 
## ... Similar to previous best
## Run 13 stress 0.2171489 
## Run 14 stress 0.1841584 
## Run 15 stress 0.1834317 
## ... Procrustes: rmse 0.00002080617  max resid 0.00007313082 
## ... Similar to previous best
## Run 16 stress 0.1834317 
## ... Procrustes: rmse 0.000004696044  max resid 0.00001909514 
## ... Similar to previous best
## Run 17 stress 0.1841584 
## Run 18 stress 0.1841584 
## Run 19 stress 0.1834317 
## ... Procrustes: rmse 0.000008098774  max resid 0.00002671289 
## ... Similar to previous best
## Run 20 stress 0.1834317 
## ... New best solution
## ... Procrustes: rmse 0.000006132046  max resid 0.00002011371 
## ... Similar to previous best
## *** Best solution repeated 1 times

The first plot is clr normalization (compositional), the second is vst(DESeq2). In contrast to bacteria - for vst WGCNA did not show good results (no plateau yield).
We still applied vst stabilization, but we have to keep in mind that the resulting output is not great.

otu.inall <- as.data.frame(ps.its.inall@otu_table@.Data) 
ps.inall.clr <- ps.its.inall
otu_table(ps.inall.clr) <- phyloseq::otu_table(compositions::clr(otu.inall), taxa_are_rows = FALSE)
data <- ps.inall.clr@otu_table@.Data %>% 
  as.data.frame()
rownames(data) <- as.character(ps.inall.clr@sam_data$Description)
powers <-  c(c(1:10), seq(from = 12, to=30, by=1))
sft <-  pickSoftThreshold(data, powerVector = powers, verbose = 5, networkType = "signed")
## pickSoftThreshold: will use block size 116.
##  pickSoftThreshold: calculating connectivity for given powers...
##    ..working on genes 1 through 116 of 116
##    Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
## 1      1   0.3860 -25.70         0.5550 57.2000 57.100000 59.600
## 2      2   0.4580 -17.60         0.3190 29.5000 29.300000 32.400
## 3      3   0.4860 -11.60         0.3420 15.8000 15.600000 18.500
## 4      4   0.6370  -7.48         0.5340  8.7100  8.530000 11.300
## 5      5   0.7900  -5.30         0.7300  4.9900  4.780000  7.370
## 6      6   0.8050  -3.89         0.7500  2.9600  2.770000  5.090
## 7      7   0.8400  -2.94         0.8030  1.8300  1.650000  3.710
## 8      8   0.7890  -2.65         0.7530  1.1700  1.000000  2.820
## 9      9   0.2080  -7.17        -0.0165  0.7790  0.626000  2.230
## 10    10   0.8000  -1.88         0.8250  0.5400  0.401000  1.810
## 11    12   0.8610  -1.51         0.9500  0.2890  0.173000  1.360
## 12    13   0.7990  -1.44         0.8150  0.2230  0.121000  1.220
## 13    14   0.8000  -1.34         0.8260  0.1770  0.083400  1.120
## 14    15   0.7480  -1.32         0.7690  0.1440  0.057500  1.020
## 15    16   0.1540  -2.22        -0.0813  0.1190  0.040100  0.941
## 16    17   0.1530  -2.10        -0.0816  0.1010  0.027800  0.869
## 17    18   0.8570  -1.15         0.8780  0.0866  0.019400  0.805
## 18    19   0.0897  -1.46        -0.0687  0.0754  0.013700  0.746
## 19    20   0.0930  -1.43        -0.0720  0.0663  0.009910  0.693
## 20    21   0.2140  -2.01         0.0509  0.0589  0.007180  0.660
## 21    22   0.2200  -2.04         0.0648  0.0527  0.005230  0.643
## 22    23   0.0951  -1.76        -0.0800  0.0475  0.003800  0.627
## 23    24   0.1930  -1.80        -0.0248  0.0431  0.002730  0.612
## 24    25   0.1960  -1.82        -0.0191  0.0393  0.001970  0.598
## 25    26   0.2010  -2.33        -0.0126  0.0360  0.001440  0.585
## 26    27   0.2220  -2.37         0.0609  0.0331  0.001070  0.572
## 27    28   0.2290  -2.40         0.0789  0.0305  0.000788  0.559
## 28    29   0.2320  -2.67         0.0387  0.0283  0.000584  0.547
## 29    30   0.2430  -2.73         0.0567  0.0263  0.000433  0.536
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n", main = paste("Scale independence"))
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], labels=powers,cex=0.9,col="red")
abline(h=0.9,col="salmon")

ps.varstab <- ps_vst(ps.its.inall, "Day")

data2 <- ps.varstab@otu_table@.Data %>% 
  as.data.frame()

rownames(data2) <- as.character(ps.varstab@sam_data$Description)
powers <-  c(seq(from = 1, to=10, by=0.5), seq(from = 11, to=20, by=1))
sft2 <-  pickSoftThreshold(data2, powerVector = powers, verbose = 5, networkType = "signed hybrid")
## pickSoftThreshold: will use block size 116.
##  pickSoftThreshold: calculating connectivity for given powers...
##    ..working on genes 1 through 116 of 116
##    Power SFT.R.sq  slope truncated.R.sq mean.k.  median.k. max.k.
## 1    1.0    0.352 -1.180        0.82400 11.7000 11.1000000  21.40
## 2    1.5    0.511 -0.965        0.80800  6.4500  5.9700000  13.20
## 3    2.0    0.688 -1.020        0.82800  3.9000  3.3000000   9.04
## 4    2.5    0.787 -0.915        0.80900  2.5200  1.9300000   6.41
## 5    3.0    0.828 -0.797        0.78200  1.7200  1.1700000   4.69
## 6    3.5    0.807 -0.913        0.75900  1.2300  0.7570000   3.81
## 7    4.0    0.789 -1.110        0.73200  0.9160  0.5010000   3.54
## 8    4.5    0.200 -2.790        0.08170  0.7060  0.3370000   3.34
## 9    5.0    0.218 -2.720        0.11300  0.5610  0.2300000   3.18
## 10   5.5    0.235 -3.520        0.08280  0.4570  0.1590000   3.04
## 11   6.0    0.165 -2.140       -0.07370  0.3800  0.1100000   2.92
## 12   6.5    0.204 -2.250       -0.00278  0.3220  0.0758000   2.82
## 13   7.0    0.208 -2.210        0.00364  0.2780  0.0538000   2.73
## 14   7.5    0.192 -2.760       -0.03940  0.2430  0.0391000   2.65
## 15   8.0    0.192 -2.670       -0.03910  0.2150  0.0280000   2.57
## 16   8.5    0.196 -2.600       -0.03270  0.1920  0.0203000   2.50
## 17   9.0    0.228 -2.770        0.01890  0.1740  0.0148000   2.44
## 18   9.5    0.213 -2.990       -0.00815  0.1580  0.0108000   2.37
## 19  10.0    0.214 -3.060       -0.00549  0.1450  0.0078000   2.32
## 20  11.0    0.230 -3.170        0.01010  0.1240  0.0041500   2.21
## 21  12.0    0.262 -3.210        0.05580  0.1080  0.0024100   2.12
## 22  13.0    0.222 -3.130        0.01010  0.0955  0.0014200   2.03
## 23  14.0    0.234 -3.060        0.02160  0.0855  0.0008090   1.95
## 24  15.0    0.259 -3.140        0.04770  0.0774  0.0004490   1.88
## 25  16.0    0.270 -3.070        0.06090  0.0707  0.0002500   1.82
## 26  17.0    0.233 -2.840        0.04610  0.0651  0.0001400   1.76
## 27  18.0    0.244 -2.790        0.05240  0.0602  0.0000789   1.71
## 28  19.0    0.251 -2.860        0.06040  0.0561  0.0000444   1.66
## 29  20.0    0.280 -3.030        0.08580  0.0525  0.0000250   1.62
plot(sft2$fitIndices[,1], -sign(sft2$fitIndices[,3])*sft2$fitIndices[,2], xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n", main = paste("Scale independence"))
text(sft2$fitIndices[,1], -sign(sft2$fitIndices[,3])*sft2$fitIndices[,2], labels=powers,cex=0.9,col="red")
abline(h=0.9,col="salmon")

# WGCNA, as old and not supported
detachAllPackages()
library(WGCNA)

net3 <- WGCNA::blockwiseModules(data2,
                          power=3,
                          TOMType="signed",
                          networkType="signed hybrid",
                          nThreads=0)


mergedColors2 <-  WGCNA::labels2colors(net3$colors, colorSeq = c("salmon", "darkgreen", "cyan", "red", "blue", "plum"))

plotDendroAndColors(
  net3$dendrograms[[1]],
  mergedColors2[net3$blockGenes[[1]]],
  "Module colors",
  dendroLabels = FALSE,
  hang = 0.03,
  addGuide = TRUE,
  guideHang = 0.05)

library(phyloseq)
library(tidyverse)
library(ggpubr)
library(ampvis2)
library(heatmaply)
library(WGCNA)
library(phyloseq)
library(ggtree)
library(tidyverse)
library(KneeArrower)

The green cluster is the background, which in the case of the red cluster for bacteria that here, I don’t think it can be called a cluster at all.
WGCNA as an analysis helps clustering based on correlations by finding a soft threshold(power, see plot where the red numbers still can’t reach the plateau),
the darkgreen and red cluster are simply unconnected phylotypes.

modules_of_interest = mergedColors2 %>% 
  unique()

module_df <- data.frame(
  asv = names(net3$colors),
  colors = mergedColors2
)

# module_df[module_df == "yellow"] <- "salmon"

submod <-  module_df %>%
  subset(colors %in% modules_of_interest)

row.names(module_df) = module_df$asv

subexpr = as.data.frame(t(data2))[submod$asv,]

submod_df <- data.frame(subexpr) %>%
  mutate(
    asv = row.names(.)
  ) %>%
  pivot_longer(-asv) %>%
  mutate(
    module = module_df[asv,]$colors
  )

submod_df <- submod_df %>% 
  mutate(name = gsub("\\_.*","",submod_df$name)) %>% 
  group_by(name, asv) %>% 
  summarise(value = mean(value), asv = asv, module = module) %>% 
  relocate(c(asv, name, value, module)) %>% 
  ungroup() %>% 
  mutate(module = as.factor(module))
## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
##   always returns an ungrouped data frame and adjust accordingly.
p <- submod_df %>% 
  ggplot(., aes(x=name, y=value, group=asv)) +
  geom_line(aes(color = module),
            alpha = 0.2) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 90),
    legend.position = "none") +
  facet_grid(rows = vars(module)) +
  labs(x = "day",
       y = "normalized abundance")

p + scale_color_manual(values = levels(submod_df$module)) +
  scale_x_discrete(labels=(c("D01"="3",
                                "D03"="14",
                                "D05"="28",
                                "D07"="49" ,
                                "D08"="63",
                                "D10"="91",
                                "D12"="119",
                                "D13"="140",
                                "D14"="161",
                                "D15"="182")
         ))

l_its <- color_filt_broken(ps_its, submod_df, ps.its.inall)
l_its

$cyan \(cyan\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 36 taxa and 36 samples ] sample_data() Sample Data: [ 36 samples by 3 sample variables ] tax_table() Taxonomy Table: [ 36 taxa by 7 taxonomic ranks ] refseq() DNAStringSet: [ 36 reference sequences ]

\(cyan\)amp ampvis2 object with 4 elements. Summary of OTU table: [1] 36 36 318660 0 23584 9007.5 8851.67

Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 28(78%) 25(69.44%) 17(47.22%) 12(33.33%) 12(33.33%) 12(33.33%) 9(25%)

Metadata variables: 4 SampleID, Day, Description, Group

\(cyan\)heat \(cyan\)heat_rel \(cyan\)tree \(cyan\)taxa

Kingdom Phylum Class Order Family Genus Species
Seq1 k__Fungi p__Ascomycota c__Sordariomycetes o__Chaetosphaeriales f__Chaetosphaeriaceae g__Chloridium s__aseptatum
Seq2 k__Fungi p__Ascomycota c__Sordariomycetes o__Coniochaetales f__Coniochaetaceae g__Coniochaeta s__verticillata
Seq12 k__Fungi p__Ascomycota c__Sordariomycetes NA NA NA NA
Seq184 NA NA NA NA NA NA NA
Seq10 k__Fungi p__Ascomycota c__Sordariomycetes o__Coniochaetales f__Coniochaetaceae g__Coniochaeta s__canina
Seq73 k__Alveolata NA NA NA NA NA NA
Seq5 k__Fungi p__Ascomycota c__Sordariomycetes o__Coniochaetales f__Coniochaetaceae g__Coniochaeta s__verticillata
Seq118 k__Metazoa p__Nematoda NA NA NA NA NA
Seq18 k__Fungi p__Ascomycota c__Sordariomycetes NA NA NA NA
Seq93 k__Fungi p__Ascomycota c__Dothideomycetes o__Venturiales f__Sympoventuriaceae g__Scolecobasidium s__constrictum
Seq78 k__Metazoa p__Nematoda NA NA NA NA NA
Seq101 NA NA NA NA NA NA NA
Seq4 k__Fungi p__Ascomycota c__Sordariomycetes o__Coniochaetales f__Coniochaetaceae g__Coniochaeta s__verticillata
Seq62 k__Metazoa p__Nematoda NA NA NA NA NA
Seq48 k__Metazoa p__Nematoda NA NA NA NA NA
Seq36 k__Fungi p__Ascomycota c__Sordariomycetes NA NA NA NA
Seq44 k__Metazoa p__Nematoda NA NA NA NA NA
Seq103 k__Metazoa p__Nematoda NA NA NA NA NA
Seq166 k__Fungi p__Ascomycota c__Eurotiomycetes o__Chaetothyriales f__Herpotrichiellaceae g__Exophiala NA
Seq231 k__Eukaryota_kgd_Incertae_sedis NA NA NA NA NA NA
Seq21 k__Fungi p__Ascomycota c__Sordariomycetes o__Hypocreales f__Clavicipitaceae g__Metarhizium s__marquandii
Seq108 k__Metazoa p__Nematoda NA NA NA NA NA
Seq74 k__Fungi p__Ascomycota c__Sordariomycetes NA NA NA NA
Seq195 NA NA NA NA NA NA NA
Seq28 k__Fungi p__Ascomycota c__Leotiomycetes o__Helotiales f__Helotiaceae g__Scytalidium NA
Seq104 k__Metazoa p__Nematoda c__Chromadorea o__Rhabditida f__Cephalobidae g__Pseudacrobeles NA
Seq29 k__Fungi p__Ascomycota c__Sordariomycetes o__Glomerellales f__Glomerellaceae g__Colletotrichum s__sidae
Seq46 k__Fungi p__Ascomycota c__Sordariomycetes NA NA NA NA
Seq122 NA NA NA NA NA NA NA
Seq235 NA NA NA NA NA NA NA
Seq191 NA NA NA NA NA NA NA
Seq94 k__Alveolata NA NA NA NA NA NA
Seq76 k__Fungi p__Ascomycota c__Sordariomycetes o__Sordariales f__Sordariales_fam_Incertae_sedis g__Staphylotrichum s__boninense
Seq155 NA NA NA NA NA NA NA
Seq190 NA NA NA NA NA NA NA
Seq32 k__Metazoa p__Nematoda NA NA NA NA NA

$darkgreen \(darkgreen\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 56 taxa and 36 samples ] sample_data() Sample Data: [ 36 samples by 3 sample variables ] tax_table() Taxonomy Table: [ 56 taxa by 7 taxonomic ranks ] refseq() DNAStringSet: [ 56 reference sequences ]

\(darkgreen\)amp ampvis2 object with 4 elements. Summary of OTU table: [1] 36 56 261584 1790 20315 6577 7266.22

Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 52(93%) 51(91.07%) 48(85.71%) 48(85.71%) 48(85.71%) 48(85.71%) 34(60.71%)

Metadata variables: 4 SampleID, Day, Description, Group

\(darkgreen\)heat \(darkgreen\)heat_rel \(darkgreen\)tree \(darkgreen\)taxa

Kingdom Phylum Class Order Family Genus Species
Seq92 k__Fungi p__Basidiomycota c__Agaricomycetes o__Agaricales f__Bolbitiaceae g__Conocybe s__zeylanica
Seq49 k__Fungi p__Basidiomycota c__Agaricomycetes o__Agaricales f__Bolbitiaceae g__Conocybe s__zeylanica
Seq91 k__Fungi p__Ascomycota c__Sordariomycetes o__Coniochaetales f__Coniochaetaceae g__Coniochaeta s__canina
Seq6 k__Fungi p__Ascomycota c__Sordariomycetes o__Sordariales f__Chaetomiaceae g__Humicola s__sardiniae
Seq16 k__Fungi p__Basidiomycota c__Agaricomycetes o__Agaricales f__Bolbitiaceae g__Conocybe s__zeylanica
Seq67 k__Metazoa p__Annelida c__Clitellata o__Enchytraeida f__Enchytraeidae g__Fridericia NA
Seq13 k__Fungi p__Ascomycota c__Eurotiomycetes o__Eurotiales f__Trichocomaceae g__Talaromyces NA
Seq7 k__Fungi p__Basidiomycota c__Agaricomycetes o__Agaricales f__Bolbitiaceae g__Conocybe s__zeylanica
Seq25 k__Fungi p__Basidiomycota c__Cystobasidiomycetes o__Cystobasidiales f__Cystobasidiaceae g__Occultifur NA
Seq69 k__Fungi p__Ascomycota c__Sordariomycetes o__Hypocreales f__Nectriaceae g__Gibberella s__intricans
Seq86 k__Fungi p__Basidiomycota c__Agaricomycetes o__Agaricales f__Bolbitiaceae g__Conocybe s__zeylanica
Seq205 k__Fungi p__Mucoromycota c__Mucoromycetes o__Mucorales f__Rhizopodaceae g__Rhizopus s__arrhizus
Seq50 k__Fungi p__Basidiomycota c__Cystobasidiomycetes o__Cystobasidiales f__Cystobasidiaceae g__Occultifur NA
Seq193 k__Heterolobosa p__Heterolobosa_phy_Incertae_sedis c__Heterolobosea o__Schizopyrenida f__Vahlkampfiidae g__Naegleria NA
Seq136 k__Fungi p__Ascomycota c__Eurotiomycetes o__Chaetothyriales f__Herpotrichiellaceae g__Exophiala NA
Seq23 k__Fungi p__Mucoromycota c__Mucoromycetes o__Mucorales f__Mucoraceae g__Actinomucor s__elegans
Seq175 k__Fungi p__Ascomycota c__Sordariomycetes o__Hypocreales f__Bionectriaceae g__Clonostachys s__rosea
Seq281 k__Heterolobosa p__Heterolobosa_phy_Incertae_sedis c__Heterolobosea o__Schizopyrenida f__Vahlkampfiidae g__Allovahlkampfia NA
Seq47 k__Fungi p__Ascomycota c__Sordariomycetes o__Hypocreales f__Nectriaceae g__Gibberella s__intricans
Seq137 k__Fungi p__Ascomycota c__Dothideomycetes o__Venturiales f__Sympoventuriaceae g__Ochroconis s__tshawytschae
Seq9 k__Fungi p__Ascomycota c__Sordariomycetes o__Hypocreales f__Stachybotryaceae g__Albifimbria s__verrucaria
Seq42 k__Fungi p__Basidiomycota c__Agaricomycetes o__Agaricales f__Bolbitiaceae g__Conocybe s__zeylanica
Seq45 k__Fungi NA NA NA NA NA NA
Seq11 k__Fungi p__Ascomycota c__Sordariomycetes o__Hypocreales f__Nectriaceae g__Gibberella s__intricans
Seq89 k__Fungi p__Basidiomycota c__Cystobasidiomycetes o__Cystobasidiales f__Cystobasidiaceae g__Occultifur NA
Seq171 k__Fungi p__Ascomycota c__Dothideomycetes o__Venturiales f__Sympoventuriaceae g__Ochroconis s__tshawytschae
Seq178 k__Fungi p__Ascomycota c__Dothideomycetes o__Venturiales f__Sympoventuriaceae g__Scolecobasidium s__constrictum
Seq22 k__Fungi p__Ascomycota c__Sordariomycetes o__Coniochaetales f__Coniochaetaceae g__Coniochaeta s__verticillata
Seq31 k__Fungi p__Mucoromycota c__Mucoromycetes o__Mucorales f__Mucoraceae g__Actinomucor s__elegans
Seq65 k__Fungi p__Ascomycota c__Sordariomycetes o__Coniochaetales f__Coniochaetaceae g__Coniochaeta NA
Seq75 NA NA NA NA NA NA NA
Seq85 k__Metazoa p__Nematoda NA NA NA NA NA
Seq54 k__Fungi p__Ascomycota c__Sordariomycetes o__Coniochaetales f__Coniochaetaceae g__Coniochaeta NA
Seq121 k__Fungi p__Ascomycota c__Sordariomycetes o__Hypocreales f__Nectriaceae g__Gibberella s__intricans
Seq125 k__Fungi p__Ascomycota c__Dothideomycetes o__Venturiales f__Sympoventuriaceae g__Ochroconis s__tshawytschae
Seq181 k__Fungi p__Ascomycota c__Eurotiomycetes o__Eurotiales f__Aspergillaceae g__Penicillium s__bialowiezense
Seq24 k__Fungi p__Ascomycota c__Sordariomycetes o__Hypocreales f__Nectriaceae g__Gibberella s__intricans
Seq27 k__Fungi p__Ascomycota c__Sordariomycetes o__Hypocreales f__Stachybotryaceae g__Stachybotrys s__chartarum
Seq66 k__Fungi p__Ascomycota c__Sordariomycetes o__Hypocreales f__Stachybotryaceae g__Albifimbria s__verrucaria
Seq114 k__Fungi p__Basidiomycota c__Agaricomycetes o__Agaricales f__Bolbitiaceae g__Conocybe s__zeylanica
Seq15 k__Fungi p__Ascomycota c__Sordariomycetes o__Hypocreales f__Nectriaceae g__Gibberella s__intricans
Seq226 k__Fungi p__Ascomycota c__Dothideomycetes o__Venturiales f__Sympoventuriaceae g__Scolecobasidium s__constrictum
Seq26 k__Fungi p__Ascomycota c__Sordariomycetes o__Sordariales f__Chaetomiaceae g__Chaetomium s__iranianum
Seq59 k__Fungi p__Ascomycota c__Dothideomycetes o__Venturiales f__Sympoventuriaceae g__Ochroconis s__tshawytschae
Seq41 k__Fungi p__Basidiomycota c__Agaricomycetes o__Cantharellales f__Ceratobasidiaceae g__Waitea s__circinata
Seq157 k__Fungi p__Ascomycota c__Sordariomycetes o__Sordariales f__Chaetomiaceae g__Zopfiella NA
Seq124 k__Fungi p__Ascomycota c__Eurotiomycetes o__Eurotiales f__Aspergillaceae g__Aspergillus NA
Seq254 k__Metazoa p__Nematoda NA NA NA NA NA
Seq17 k__Metazoa p__Nematoda NA NA NA NA NA
Seq39 k__Fungi p__Ascomycota c__Sordariomycetes o__Coniochaetales f__Coniochaetaceae g__Coniochaeta s__verticillata
Seq180 NA NA NA NA NA NA NA
Seq20 k__Fungi p__Ascomycota c__Sordariomycetes o__Sordariales f__Chaetomiaceae g__Chaetomium NA
Seq102 NA NA NA NA NA NA NA
Seq70 k__Fungi p__Ascomycota c__Sordariomycetes o__Sordariales f__Sordariales_fam_Incertae_sedis g__Staphylotrichum s__boninense
Seq150 k__Heterolobosa p__Heterolobosa_phy_Incertae_sedis c__Heterolobosea o__Schizopyrenida f__Vahlkampfiidae g__Allovahlkampfia NA
Seq132 NA NA NA NA NA NA NA

$salmon \(salmon\)ps phyloseq-class experiment-level object otu_table() OTU Table: [ 24 taxa and 36 samples ] sample_data() Sample Data: [ 36 samples by 3 sample variables ] tax_table() Taxonomy Table: [ 24 taxa by 7 taxonomic ranks ] refseq() DNAStringSet: [ 24 reference sequences ]

\(salmon\)amp ampvis2 object with 4 elements. Summary of OTU table: [1] 36 24 121335 20 14955 1772 3370.42

Assigned taxonomy: Kingdom Phylum Class Order Family Genus Species 22(92%) 18(75%) 15(62.5%) 14(58.33%) 13(54.17%) 13(54.17%) 11(45.83%)

Metadata variables: 4 SampleID, Day, Description, Group

\(salmon\)heat \(salmon\)heat_rel \(salmon\)tree \(salmon\)taxa

Kingdom Phylum Class Order Family Genus Species
Seq3 k__Fungi p__Ascomycota c__Sordariomycetes o__Sordariales f__Lasiosphaeriaceae g__Schizothecium s__inaequale
Seq19 k__Fungi p__Basidiomycota c__Agaricomycetes o__Agaricales f__Psathyrellaceae g__Coprinellus s__flocculosus
Seq53 k__Fungi p__Ascomycota c__Sordariomycetes o__Hypocreales f__Stachybotryaceae g__Albifimbria s__verrucaria
Seq99 k__Viridiplantae p__Anthophyta NA NA NA NA NA
Seq215 k__Viridiplantae NA NA NA NA NA NA
Seq30 k__Fungi p__Ascomycota c__Sordariomycetes o__Sordariales f__Chaetomiaceae g__Humicola s__sardiniae
Seq43 k__Fungi p__Ascomycota c__Sordariomycetes o__Sordariales f__Lasiosphaeriaceae g__Schizothecium NA
Seq68 NA NA NA NA NA NA NA
Seq63 k__Fungi p__Mucoromycota c__Mucoromycetes o__Mucorales f__Mucoraceae g__Actinomucor s__elegans
Seq100 k__Fungi NA NA NA NA NA NA
Seq140 k__Fungi p__Ascomycota c__Dothideomycetes o__Pleosporales NA NA NA
Seq33 k__Fungi p__Basidiomycota c__Agaricomycetes o__Agaricales f__Psathyrellaceae g__Coprinellus s__flocculosus
Seq245 NA NA NA NA NA NA NA
Seq138 k__Viridiplantae p__Anthophyta c__Eudicotyledonae NA NA NA NA
Seq8 k__Fungi p__Ascomycota c__Sordariomycetes o__Sordariales f__Lasiosphaeriaceae g__Schizothecium s__inaequale
Seq79 k__Viridiplantae p__Anthophyta NA NA NA NA NA
Seq161 k__Eukaryota_kgd_Incertae_sedis NA NA NA NA NA NA
Seq96 k__Metazoa p__Nematoda c__Chromadorea o__Rhabditida f__Cephalobidae g__Acrobeloides s__nanus
Seq147 k__Fungi NA NA NA NA NA NA
Seq123 k__Fungi p__Ascomycota c__Sordariomycetes o__Sordariales f__Chaetomiaceae g__Chaetomium s__jodhpurense
Seq152 k__Viridiplantae p__Anthophyta NA NA NA NA NA
Seq52 k__Fungi p__Ascomycota c__Sordariomycetes o__Sordariales f__Lasiosphaeriaceae g__Schizothecium NA
Seq40 k__Fungi p__Mucoromycota c__Mucoromycetes o__Mucorales f__Mucoraceae g__Actinomucor s__elegans
Seq37 k__Fungi p__Ascomycota c__Sordariomycetes o__Hypocreales f__Nectriaceae g__Gibberella s__intricans
p.observed <- plot_alpha_w_toc(ps = ps_its, group = "Group", metric = c("Observed")) + 
  theme(axis.title.y = element_blank())

p.shannon <- plot_alpha_w_toc(ps = ps_its, group = "Group", metric = c("Shannon")) + 
  theme(axis.title.y = element_blank())

p.simpson <- plot_alpha_w_toc(ps = ps_its, group = "Group", metric = c("InvSimpson")) + 
  theme(axis.title.y = element_blank())

ggpubr::ggarrange(p.observed, p.shannon, p.simpson, ncol = 3)

p.observed <- plot_alpha_w_toc(ps = ps_its, group = "Day", metric = c("Observed")) + 
  theme(axis.title.y = element_blank())

p.shannon <- plot_alpha_w_toc(ps = ps_its, group = "Day", metric = c("Shannon")) + 
  theme(axis.title.y = element_blank())

p.simpson <- plot_alpha_w_toc(ps = ps_its, group = "Day", metric = c("InvSimpson")) + 
  theme(axis.title.y = element_blank())

ggpubr::ggarrange(p.observed, p.shannon, p.simpson, ncol = 3)

ps.m <- phyloseq::psmelt(ps_its)
ps.m <- ps.m %>% 
  mutate_if(is.character, as.factor) 

ps.data.out <- ps.m %>%
  select(-Group) %>% 
  pivot_wider(names_from = c(Day, Description, Sample), values_from = Abundance, values_fill = 0) 

#create empty dataframe with columnnames
external_empty_dataframe <- data.frame(OTU=factor(), cluster=factor(), stringsAsFactors = FALSE)

for (i in names(l_its)) {
  a <- taxa_names(l_its[[i]][["ps"]])
  b <- rep(i, length(a))
  d <- data.frame(OTU = as.factor(a),
                  cluster = as.factor(b))
  external_empty_dataframe <- rbind(external_empty_dataframe, d)
}

clusters.otu.df <- external_empty_dataframe
# add exl taxa -- taxa "exl"
d <- data.frame(OTU = as.factor(ps.its.exl.taxa),
                cluster = as.factor(rep("exl", length(ps.its.exl.taxa))))

clusters.otu.df <- rbind(clusters.otu.df, d)
ps.data.out.exl <- left_join(clusters.otu.df, ps.data.out, by="OTU")
# write.table(ps.data.out.exl, file = "ps.its.data.out.tsv", sep = "\t")

0.12 Packages info

sessionInfo()
## R version 4.2.0 (2022-04-22)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.6 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] gridExtra_2.3         reshape2_1.4.4        KneeArrower_1.0.0    
##  [4] ggtree_3.6.2          heatmaply_1.4.2       viridis_0.6.2        
##  [7] viridisLite_0.4.1     plotly_4.10.1         ampvis2_2.7.28       
## [10] ggpubr_0.6.0          lubridate_1.9.2       forcats_1.0.0        
## [13] stringr_1.5.0         dplyr_1.1.0           purrr_1.0.1          
## [16] readr_2.1.4           tidyr_1.3.0           tibble_3.1.8         
## [19] ggplot2_3.4.1         tidyverse_2.0.0       phyloseq_1.42.0      
## [22] WGCNA_1.72-1          fastcluster_1.2.3     dynamicTreeCut_1.63-1
## 
## loaded via a namespace (and not attached):
##   [1] estimability_1.4.1             coda_0.19-4                   
##   [3] bit64_4.0.5                    knitr_1.42                    
##   [5] irlba_2.3.5.1                  multcomp_1.4-22               
##   [7] DelayedArray_0.24.0            data.table_1.14.8             
##   [9] rpart_4.1.19                   KEGGREST_1.38.0               
##  [11] RCurl_1.98-1.10                doParallel_1.0.17             
##  [13] generics_0.1.3                 preprocessCore_1.60.2         
##  [15] BiocGenerics_0.44.0            ScaledMatrix_1.6.0            
##  [17] cowplot_1.1.1                  TH.data_1.1-1                 
##  [19] RSQLite_2.3.0                  proxy_0.4-27                  
##  [21] bit_4.0.5                      tzdb_0.3.0                    
##  [23] bayesm_3.1-5                   webshot_0.5.4                 
##  [25] SummarizedExperiment_1.28.0    assertthat_0.2.1              
##  [27] DirichletMultinomial_1.40.0    xfun_0.37                     
##  [29] hms_1.1.2                      jquerylib_0.1.4               
##  [31] evaluate_0.20                  TSP_1.2-2                     
##  [33] DEoptimR_1.0-11                fansi_1.0.4                   
##  [35] dendextend_1.16.0              readxl_1.4.2                  
##  [37] geneplotter_1.76.0             mia_1.6.0                     
##  [39] igraph_1.4.1                   DBI_1.1.3                     
##  [41] htmlwidgets_1.6.1              tensorA_0.36.2                
##  [43] stats4_4.2.0                   Rmpfr_0.9-1                   
##  [45] CVXR_1.0-11                    ellipsis_0.3.2                
##  [47] corrplot_0.92                  crosstalk_1.2.0               
##  [49] backports_1.4.1                energy_1.7-11                 
##  [51] signal_0.7-7                   permute_0.9-7                 
##  [53] picante_1.8.2                  annotate_1.76.0               
##  [55] compositions_2.0-5             deldir_1.0-6                  
##  [57] sparseMatrixStats_1.10.0       MatrixGenerics_1.10.0         
##  [59] vctrs_0.5.2                    SingleCellExperiment_1.20.0   
##  [61] Biobase_2.58.0                 abind_1.4-5                   
##  [63] cachem_1.0.7                   withr_2.5.0                   
##  [65] ggforce_0.4.1                  robustbase_0.95-0             
##  [67] vroom_1.6.1                    checkmate_2.1.0               
##  [69] emmeans_1.8.4-1                vegan_2.6-4                   
##  [71] treeio_1.22.0                  MultiAssayExperiment_1.24.0   
##  [73] mnormt_2.1.1                   cluster_2.1.4                 
##  [75] gsl_2.1-8                      ape_5.7                       
##  [77] lazyeval_0.2.2                 crayon_1.5.2                  
##  [79] TreeSummarizedExperiment_2.6.0 labeling_0.4.2                
##  [81] pkgconfig_2.0.3                tweenr_2.0.2                  
##  [83] GenomeInfoDb_1.34.9            nlme_3.1-162                  
##  [85] vipor_0.4.5                    seriation_1.4.1               
##  [87] nnet_7.3-18                    rlang_1.0.6                   
##  [89] lifecycle_1.0.3                sandwich_3.0-2                
##  [91] registry_0.5-1                 rsvd_1.0.5                    
##  [93] polyclip_1.10-4                cellranger_1.1.0              
##  [95] matrixStats_0.63.0             rngtools_1.5.2                
##  [97] aplot_0.1.9                    Matrix_1.5-3                  
##  [99] carData_3.0-5                  Rhdf5lib_1.20.0               
## [101] boot_1.3-28.1                  zoo_1.8-11                    
## [103] base64enc_0.1-3                beeswarm_0.4.0                
## [105] png_0.1-8                      ca_0.71.1                     
## [107] rootSolve_1.8.2.3              bitops_1.0-7                  
## [109] rhdf5filters_1.10.0            Biostrings_2.66.0             
## [111] blob_1.2.3                     DelayedMatrixStats_1.20.0     
## [113] doRNG_1.8.6                    decontam_1.18.0               
## [115] gridGraphics_0.5-1             jpeg_0.1-10                   
## [117] rstatix_0.7.2                  DECIPHER_2.26.0               
## [119] S4Vectors_0.36.2               ggsignif_0.6.4                
## [121] beachmat_2.14.0                scales_1.2.1                  
## [123] memoise_2.0.1                  magrittr_2.0.3                
## [125] plyr_1.8.8                     zlibbioc_1.44.0               
## [127] compiler_4.2.0                 RColorBrewer_1.1-3            
## [129] lme4_1.1-31                    DESeq2_1.38.3                 
## [131] cli_3.6.0                      ade4_1.7-22                   
## [133] ANCOMBC_2.0.2                  XVector_0.38.0                
## [135] lmerTest_3.1-3                 patchwork_1.1.2               
## [137] htmlTable_2.4.1                Formula_1.2-5                 
## [139] MASS_7.3-58.2                  mgcv_1.8-42                   
## [141] tidyselect_1.2.0               stringi_1.7.12                
## [143] highr_0.10                     yaml_2.3.7                    
## [145] locfit_1.5-9.7                 BiocSingular_1.14.0           
## [147] latticeExtra_0.6-30            ggrepel_0.9.3                 
## [149] grid_4.2.0                     sass_0.4.5                    
## [151] tools_4.2.0                    lmom_2.9                      
## [153] timechange_0.2.0               parallel_4.2.0                
## [155] rstudioapi_0.14                foreach_1.5.2                 
## [157] foreign_0.8-84                 gld_2.6.6                     
## [159] farver_2.1.1                   digest_0.6.31                 
## [161] Rcpp_1.0.10                    GenomicRanges_1.50.2          
## [163] car_3.1-1                      broom_1.0.3                   
## [165] scuttle_1.8.4                  AnnotationDbi_1.60.0          
## [167] httr_1.4.5                     psych_2.2.9                   
## [169] Rdpack_2.4                     colorspace_2.1-0              
## [171] XML_3.99-0.13                  IRanges_2.32.0                
## [173] splines_4.2.0                  yulab.utils_0.0.6             
## [175] tidytree_0.4.2                 expm_0.999-7                  
## [177] scater_1.26.1                  multtest_2.54.0               
## [179] Exact_3.2                      ggplotify_0.1.0               
## [181] systemfonts_1.0.4              xtable_1.8-4                  
## [183] gmp_0.7-1                      jsonlite_1.8.4                
## [185] nloptr_2.0.3                   ggfun_0.0.9                   
## [187] R6_2.5.1                       Hmisc_4.8-0                   
## [189] pillar_1.8.1                   htmltools_0.5.4               
## [191] DT_0.27                        glue_1.6.2                    
## [193] fastmap_1.1.1                  minqa_1.2.5                   
## [195] BiocParallel_1.32.5            BiocNeighbors_1.16.0          
## [197] class_7.3-21                   codetools_0.2-19              
## [199] mvtnorm_1.1-3                  utf8_1.2.3                    
## [201] lattice_0.20-45                bslib_0.4.2                   
## [203] numDeriv_2016.8-1.1            ggbeeswarm_0.7.1              
## [205] DescTools_0.99.48              GO.db_3.16.0                  
## [207] interp_1.1-3                   survival_3.5-3                
## [209] rmarkdown_2.20                 biomformat_1.26.0             
## [211] munsell_0.5.0                  e1071_1.7-13                  
## [213] rhdf5_2.42.0                   GenomeInfoDbData_1.2.9        
## [215] iterators_1.0.14               impute_1.72.3                 
## [217] gtable_0.3.1                   rbibutils_2.2.13